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ABSTRACT 


This thesis presents an introduction to Hidden Markov models (HMM) and their 
applications to classification problems. HMMs have been used extensively to model the 
temporal structure and variability of speech and other signals in the last decade. We 
selected to write our own HMM implementation in MATLAB. We tested our software 
on a limited isolated 4-word recognition. We also applied our implementation to the 
recognition of mine-like objects buried in shallow sand, using seismo-acoustic data 
obtained from an on-going project at the Naval Postgraduate School. Initial results 
indicate that the HMM-based classifier can recognize the type of mine-like object, 
independent of the object weight with a 97% accuracy. Results also indicate that it can 
recognize the object type at different distances with a 100% accuracy. However, the 
experiments were conducted with very few data, and further work needs to be done to 
confirm these initial findings by using a larger data set. Finally, we benchmarked our 
results against those obtained using a back-propagation neural network implementation, 


which were found to be similar, but slower than the HMM-based implementation. 
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ik. INTRODUCTION 


This thesis presents an introduction to Hidden Markov models (HMM) and their 
applications to classification problems. HMMs have been used extensively to model the 
temporal structure and variability of speech and other signals in the last decade. For 
example, they have become a major player in the speech area [4, 5, 9]. We selected to 
write our own HMM implementation even though sophisticated HMM software is readily 
available on the market to better understand the basic concepts behind the theory. Our 
implementation uses MATLAB because it is a high-level language easy to work with. 
As a result, the software may not be as fast as that obtained with a compiled 
implementation but it is easy to understand, which was one of the main goals of the 
research. We tested our software on a limited isolated 4-word recognition, and we also 
applied our implementation to the recognition of mine-like objects buried in shallow 
sand, using seismo-acoustic data obtained from an on-going project headed by Prof. Muir 
from the Physics Department at the Naval Postgraduate School. Finally, we 
benchmarked our results against those obtained using a back-propagation neural network 
implementation. 

Chapter IJ introduces the concepts of HMMs in an “engineering-oriented” simple 
fashion. We illustrate the theoretical concepts with basic examples to facilitate the 
understanding of this difficult topic. We cover the three specific problems which HMM 
address, their solutions, emphasizing the computational savings obtained with the 


forward, backward, Viterbi and Baum-Welch algorithms. 


Chapter If presents the application of HMMs to a simple speech classification 
problem, using four isolated three-syllable words: “Microsoft,” “Statistics,” “Instructor” 
and “Professor.” Our goal is to show how a generic classifier can be set-up through the 
simple speech recognition example. We introduce the concept of vector quantization 
(VQ) applied to generate discrete symbols from the speech feature vectors created with 
Linear Prediction Coefficients (LPC) and energy coefficients. Two different 
implementations of VQ are considered 1) a Neural Network (NN)-based VQ scheme; and 
2) a K-means algorithm [9]. In addition, we point out the potential numerical difficulties 
which can be encountered while setting up and implementing the HMM software. 

Chapter IV considers the application of the HMM-based classifier to the 
classification of two mine-like objects buried in sand; a cylinder and a powder keg with 
weights ranging from 71kg to 290kg. Results show that recognition performance 1s good 
under various conditions of weight and distance of the object. 

Chapter V presents a back-propagation neural network classifier designed to 
recognize the two mine-like objects discussed in Chapter IV. This implementation was 
conducted on the same data to compare the performance of the two schemes. Results 
show the performances to be similar (around 95%), but the HMM-based implementation 


procedure is faster than the NN-based implementation. 


Il. INTRODUCTION TO HIDDEN MARKOV MODELS 


This chapter introduces the basic theory of Hidden Markov models (HMM), 
which are a very powerful technique for modeling the temporal structure and variability 
in speech and other applications. The understanding of HMM is very important in order 


to use them properly and achieve the best results for signal classification. 


A. HMM BACKGROUND 


Hidden Markov Model theory was first introduced in the late 1960s and early 
1970s by Baker at Carnegie-Mellon University and Jeninek and colleagues at IBM for 
speech recognition [6]. Since then they have been used extensively for speech 
applications, and also successfully to other tasks such as human face identification, lip 
and speech-reading, optical character recognition and time DNA modeling ([6] and 
references therein). The reason for this wide range of applications is the nch 
mathematical structure the HMMs are built on, yielding optimal results if used properly. 
A detailed overview of HMM theory was presented by Rabiner in the late 1980s [3, 4]. 

There are three main reasons why we may need to model a signal. First, we apply 
signal modeling to mathematically describe the signal, so that we’ll be able to process it, 
for example to denoise a speech signal. Models are also important because they let us 
describe the signal source, which doesn’t have to be directly available to the user. For 
example we cannot produce real seismic waves without an earthquake, but we can create 
models of such signals and then process them. Finally, the most important reason for the 
widely spread use of signal models is that they perform well in practice [3]. 

There are two types of signal models; deterministic and statistical. Deterministic 


modeling 1s applied when dealing with signals with known physical characteristics. For 


example a sinewave is completely specified by its frequency, phase and amplitude. 
Stochastic modeling is applied when one tries to characterize only the statistical 
properties of the signal. For example, statistical modeling may include Gaussian Poisson 
and Markov process to describe events. Usually, real applications use both deterministic 
and stochastic modeling. In this work we focus on statistical signal modeling, using the 


Hidden Markov Model (HMM). 


B. INTRODUCTION 


The HMM theory is based on the Markov Chain. We can define the Markov 
Chain as a probabilistic description of transitions between a system’s states. A state can 
be a property, or generally a condition, that a system/model might have at a particular 
instance. The HMM consists of an underlyining Markov chain describing the 
probabilistic status between the states, as that shown in Figure 2.1, which illustrates a 
three state left-right model. For example, suppose we want to model a speech signal. 
First, the signal is split into T time frames. Then, a set of parameters (such as LPC 
coefficients, energy, etc...) 1s extracted together with a set of symbols for each time 
frame. The sets of symbols represent each frame characteristics, and are called 
observations. As a result, the entire model is a sequence of symbols, and each symbol is 
a system model depicting each segment. In most cases we choose the segment length 
empirically, but sometimes adjust it so that it is large enough to contain all the 
information (usually spectral) that makes it unique, by comparison with the other 
segments, due to possible change in the signal behavior. Generally, we can assume that 
we reach an optimal number of frames when, by decreasing the frame length, we generate 


models identical to those already generated. At this point, HMMs take advantage of the 


properties existing between adjacent segments properties by addressing the following 
three problems [3]: 1) how to identify the characteristic frames; 2) how to characterize 


the relation between all successive segments; and 3) what types of properties should be 


extracted to model each segment. 


C: DISCRETE MARKOV PROCESS 


An accurate definition for the HMM according to Rabiner [3] is that “A HMM is 
a doubly stochastic process that is not observable (it is hidden), but can only be observed 
through another set of stochastic processes that produce the sequence of observed 
symbols.” 

Consider a system which exists at time t in one of N possible potential states, as 
illustrated in Figures 2.1 and 2.2, where N=3. Each of the three circles represents a state 
of the model. At a specific discrete time instant t within frame k, the model is always at 
one of those three states, and we observe the sequence Ox. Generally, given that the 
system has N states S = {s1, S2, . . . sn}, where gj is the ie state, at every time t=k, the 
model passes through the sequences of states Q = {q, q2,.. . .dr}, where q; is the state at 
time k. The model that describes al] this information is called a Markov chain. As we 
can see from the example in Figure 2.1, 

SK Ge = Agar» (2.1) 
which means that if the model at time t=k is at state, it can remain at the same state at 
time t=k+1. Note that no backward transitions are allowed in Fig 2.1. As a result, this 
model is called a “left-to-nght” model (the model is called ergotic when backward 


transitions are allowed, as illustrated in Fig 2.2). 


We define the probability aj as the transition probability from state i to state j. 
For example, a3 is the probability of going from the 2° to the 3" state, and is defined as: 
aj = Plq = jlq -1 =i], 1Sij<N, (2.2) 


with the following constrains: 


aij > 0 Wiba je 

N 

S' ay =1 vie Ca) 
jel 


since all probabilities are positive numbers and the summation of all transition 


probabilities is 1. 


aj3 


Figure 2.1 A Markov Process (Chain), left-to-rght model 


aj] a22 





a33 


Figure 2.2 A Markov Process (Chain), ergotic model 


For the discrete-time case, a system governed by known or predictable dynamics 

can be modeled using a Markov chain. The probability of observing the observation O,, 
given the model A at a time t, is determined by the state at the time preceding states [6]: 

P(OdA) = Plae = 5j,qr-1= 8, qr-2=5k,... |. (2.4) 


For a first-order Markov chain, eq (2.4) can be simplified to: 
PO ltalselsenbeelss|h lean (25) 
k=l 


Finally, we define the initial conditions, the probabilities of starting (t=1) at the i" state s; 
as: 
m=Plg,=s,) i=1,2,---N. (2.6) 
i Example 1 
Let’s evaluate the probability of the observation O = {sj, Ss), 83} for the example 


shown in Figure 2.2. So, applying eqs (2.4) and (2.5) to the example shown in Figure 2.1 


leads to: 
Pr(O|)= Plfs: si, s3/|A] 
Pr(O|A)= Pls4s2|P[s-)s.JP[s:] 
Pr(O|1) = z,4,,4,>. 
Ze Example 2: Browser Tracking. 


Next, let us consider the following example where a user wishes to track the 
browser type used by visitors to the Web. Assume: 1) only three types of browsers are 
available: Microsoft Explorer (MIE), Netscape (NET) and America Online’ s browser 


(AOL); 2) the specific browser type version 1s not taken into account. 


Table 2.1 Browser Tracking 


Browser MS Explorer5.O | AOL4.0 Netscape3.0 Netscape6. 1 
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Browser 2 | 3 | 3 | 
Observation 
Symbol 


An observation symbol is assigned to each visitor’s browser types, as illustrated in Table 


ce 






2.1. Note that we assign the same symbol (3) to both versions of Netscape, as the 
specific browser version in not tracked by the user. In addition we assume transition 
probabilities which represent probabilities of going from one state to another, i.e., 
browser, to another are available from previous statistical studies: 


0.3 0.4 0.3 
a=|0.2 0.5 0.3}. 
0.8 0.1 O.1 


Thus, the resulting observation sequence is O={1, 2, 3, 3}. Finally, let’s assume that the 
probability the first user uses the explorer browser is 71)=0.5. 


Therefore: 


P(O|A)= Pll, 2, 3, 3] A)= PLiPl2 |1]P (3 |2]P[3|3]= z,0,,45.45, 
= P(O|A)=0.5-0.2-0.1-0.1=0.001. 


3. Example 3: Urn-ball Example 

Note that the states are direct!y observable in the previous examples. However, in 
most real world cases, the state sequence that produces a given sequence of patterns 
cannot be determinated, and the model is said to be “hidden”. The most common 


example for this type of description is the urn-ball model, according to Fig 2.3: 





ai] ao? a 


Figure 2.3. The classic Urmn-ball example. There are three (N=3) urns that contain a large 
number of color balls. The colors of the balls are red, green, and blue. The boy we assigned to 
paint the balls runs out of blue paint sometime while painting and uses sky blue afterwards. Thus, 
eventhough there are 4 colors, we assign the same symbol for blue and sky blue and the number of 
symbols is three (M=3). Each time a person randomly selects an urn, he selects a ball, and 
announces the color. The process is repeated 4 (T=4) times. Note that: 1) we don’t know from 
which urn the ball came from, as we only have a color sequence, and we map those colors with the 
symbols according from Table, 2.2; 2) we assume that it is a left-to-right model for simplicity. 


Table 2.2. Urn-ball example. T, Q, C, O values parameters 


ela Laat 


Table 2.2 shows that we have T=4 observations of four balls. The number of urns 












is N=3 and each ur represents one of the 3 hidden states, because we don’t know from 
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the specific urn number each ball comes from. Further, we consider that the sky blue 
color belongs to the same class as blue, so we assign the same observation symbol (O,=2) 
for both colors, and the total number of symbols is M=3. The possible observation 


symbols are V={1, 2, 3}, and the states are Q = {q, q2, q3}. We also assign values to the 


transition probability matrix A={ajj}, and initial conditions vector I]={7;}: 


0.5 03 02 
A={aj}=| 0 0.8 0.2], Tl = {z}= {0.9, 0.1, 0}. 
peo. i 


Note that the specific upper triangular structure of A indicates the model considered here 
is left-to-right, as we cannot go backward. In addition, we assume that the probability of 
selecting the first urn is 90%. Finally, we also need the observation symbol probability 
distribution at every state, to describe the model completely. This information is 
presented in the matrix B, which contains the probabilities of being at the state j and 
observing the symbol k, such as: 

B={bk}, 1S j<SN, 1Sk<M, (2.7) 


with the following properties (like A): 


bj =O tee 

M 2.8 
S) bik =] Wj. Ce 
lea 


Note that B is a NxM matrix where N represents the number of hidden states and M the 
number of symbols. Matrix B isn’t restricted to be square. In our example, a possible B 


matrix may be defined as: 
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0.3 0.5 0.2 
B=bx=|04 04 0.2]. 
0.6 0.3 0.1 
The statistical model is completely described by the set of matrices A, B, and the 
vector m= and usually denoted as: 
}i =) Ae ae (2.9) 
Note that all rows of matrices A, B and m shown sum-up to 1, thereby providing an easy 
check on the validity of the model. Extending our example to a general case, Table 2.3 
lists the definition of all HMM parameters, as given by Rabiner [3], for a generic HMM 
model: 


Table 2.3 HMMs elements 


T: Length of the Observation Sequence (total number of clock times) 
N: Number of States in the Model 

M: Number of observation symbols 

Q={q1, q2,---, Gn}: States 

V={v, V2, ---5 Vm}: Discrete Set of Possible Symbol Observations 


A = {aij}, aij = Pr{q; @ t+1|q; @ t}: State Transition Probability 
Distribution 


B = {b;(k)}, bj(k) = Pr(O;=vi|qi=s;): Observation Symbol Probability 
Distribution in State j 


m1 ={7,}, 7, = Pr(qj @ t = 1): Initial State Distribution 





Pe 


D. THE THREE BASIC PROBLEMS FOR HMM’S 


—— 


Three basic problems of interest must be solved in order to specify the model A = 
{A, B, 7}, and use it in classification applications [3]: 

Problem 1: How to compute the probability of the observation sequence Pr(O|A), 
for O = {O;, O2,..., Or}. 

Problem 2: How to compute the most optimal state sequence I = {i), 12, .. . ,it}, 
for a given observation sequence O = {Q;, O2,..., Or} and model 2. 

Problem 3: | How to adjust the model parameters A = {A, B, 71} to maximize 
Pr(OJA). 

Problem 1 is an evaluation problem, because we want to evaluate the probability 
Pr(O|A), for the specific model. The hidden part of the model is the state sequence I that 
we attempt to find in Problem 2. Note that there are many possible state sequences, but 
only one is optimal. Finally in Problem 3, we adjust the parameters A, B, and a of the 
model A, so that the probability Pr(O|A) is maximum, i.e., we attempt to optimize the 


model A. 


E. SOLUTIONS TO THE THREE HMM PROBLEMS 
1. Problem 1: Evaluation of Pr(O|A) 


Problem 1 deals with evaluating the probability of the observation sequence O, 
given the model A, i.e. Pr(O|A). Using basic probability principles, it is the summation of 
the conditional probability Pr(O,I|A) over all possible state sequences I [3]: 


Pr(O |A)= ¥Pr(O|7,A)Pr(7| A), (2.10) 


all J 


Ks 


This evaluation first requires the definition of Pr(O,]|A), the joint probability that the 


observation O={Q), O2,...,Or} and the state sequence I={1), in, ... , ir} occur at the 


same time, given the model A [6]. Using the Bayes’ rule it is computed as: 


Pr(O,1|A)=Pr(O|1,A)Pr( | A). 
Since we assume independence of the observation, we can write: 
rE 
Pr(O|7,2)=| | Pr(O: |i, A) =bi, (O1)bi, (O2)...b:,(Or), 
t=l 
and 
Pr(J | a qt; 493i, 4i, 213 - °Q; ir 


Therefore, replacing eqs (2.10 ), and (2.13) into (2.11 ) leads to: 


Pr(O Ve S)7,bi,(O 1)ai, i bi (O2)u, I bi (O 3)---a ar. tela (ON 


2a 3 
alli 


(2.11) 


(219) 


(2.13) 


(2.14) 


As an example of this computation, Let’ s apply the above result to the urn-ball example 


considered earlier in section C and in Figure 2.3. The components of the model 


(A, B, 7) are: 
0503.02 0.3 05 0.2 
AS|0 08 0.2)\R= 04014) 02) \wee=10.9,0 Wow 
0 O 1 06 03 0.1 


The observation sequence 1s: 


OR iN. 2, 2 55 
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with the number of states N=3, clock times T=4, and the symbol dimension M=3. Recall 
that we don’t know the state scatiens™ 1.e., we cannot identify which specific urn, the 
color balls came from. Using eq (2.14) we can find the hidden state sequence I, by 
picking that which leads to the highest probability Pr(O,I|A). For example, possible state 
sequence that we can use are I={1,1,1,1} (which means that all balls come from the first 
um), or I={1,1,1,2}. Note that this model is a left-to-right model, as the state sequence 


doesn’t go backwards, thus, for instance the state sequence I={1,3,2,2} is not valid as we 


can’t go from the 3rd to the 2nd state. Therefore: 


for 1={1,1,1, 2): (first 3 balls from 1" urn, 
last one from 2™ urn) 
mirbi(O 1)aini2bix(O 2)airisbis(O3)- = ir - irDir(Or) 
= mibi(1 Jaisixbia{2 )airisbis(2)- ‘+ Qir - irbir(3) 
= 10,4101), 9,.4,2b43 
=09-0.3-05-05-05-05-03-0.2 
= 0.0010125 


for I ={3, 3, 3,3}: (all balls from 3 urn) 
tibi(O1)aisixbi2(O2)airisbis(O 3)+++@ir - urbir(Or) 
= midi] \ainirdiz(2 )airirbis(2). ** ir - sirbir(3) 
= 1D3,033b3)433b3)33b3; 
=0.9-0.6-1-0.3-1-0.3-1-01 
= 0.00486 


As we can see, the number of combinations of all possible sequences is very 
large, even though the number of the observations T is relatively small. In this problem 
(T=4), the direct computation of Pr(O|A) requires (2T-1)N’ multiplications and N'-1 


additions, which is too expensive for real applications. For instance, the number of 


SS 


computations is (2*100-1)3'°°+3'"-1 when T=100 and N=3. Therefore, the more 


efficient forward-backward procedure was introduced to solve this problem [3, 6]. 


a) Forward Procedure: 


We define the forward variable, a,(i), as the probability of the partial 
observation sequence {QO}, Oo, ... ,O,}, until time t and state q; at time t, given the model 
dX, such as: 


a.,(i) = Pr(O,,0,,....0,,i, =qA). (2.15) 


The forward algorithm defined below is used to evaluate all possible o,4(i) variables: 


Initialization: 
(i) = =,b(O, ), 1<i<N (2.16) 
Recursion: 
N 
AG) =e aannO)  yort =, 2.0 tebe) < NV (2.17) 
i=] 
Termination: 
Pr(O|A)= 5" a, (i). (2.18) 


Eq (2.18) is also know as the Baum-Welch probability. As we can see 
from the recursion, the forward procedure gets its name from the fact that a,4, is evaluated 
using the previous value of the forward variable &%. Computing Pr(O/A) using this 
procedure requires only N°T calculations instead of the 2TN'-1 required by the direct 
definition [3]. For example direct computation requires 2TN'=7.4-10'° computations, 


while the forward procedure requires N’T=1920 for N=8 and T=30. 
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a : Current position (State, Time). Total number of positions=NeT 


- Possible movement 


} : Most likely path 


Figure 2.4a Trellis Diagram for T=4, and N=3, ergotic model 
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= : Current position (State, Time). Total number of positions=NeT 


A : Possible movement (for the left-right model they are less) 


] :Most likely path 


Figure2.4b Trellis Diagram for T=4, and N=3, left-to-right 
model as used for example in Figure 2.3. Fewer possible paths 
than in ergotic model are allowed. 


The computational savings obtained in the forward procedure are 


illustrated in Figure 2.4a and 2.4b which present all possible combinations of state 
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sequences from clock time t, to ty in a trellis diagram. Figure 2.4a presents the trellis 
diagram for an ergotic model, and T=4, N=3. Each calculation of o(1) with the forward 
procedure only requires the computation of the values %.:(j), for LSj<N. Note that the 
procedure discards the routes less likely to occur, so that the probabilities through the 


previously discarded routes do not get reevaluated at the next iteration time. 


Finally, the probability Pr(O|A) is obtained by the summation in eq (2.18). 
b) Backward Procedure: 


Similarly, we define the backward variable B,(i) as the probability of the 
partial observation sequence from t+1 to the end, given q; at time t and the model A [3], 
where: 


BG EHOW, Oe ame ae (2.19) 


B.@) is computed recursively from t=T down to t=1 as follows: 


Initialization: 
B,(i) =1, l<is<N (2.20) 
Recursion: 
N 
BO=Y2,6(0,,)8..0). t=Tl,T-2,...,1 ,15isN.(2.21) 
j=l 


Note that the backward variable is issued to re-estimate the model A (Problem 3), and that 
the number of computations is decreased to N’T. 
Therefore, using the definitions of the forward and backward variables, 


and according to [6]: 


Pr(O|i)= S a,(i)B,(i) f=. (2.22) 
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Um-Ball Example: 


Forward and backward variables are evaluated for the urn-ball example 


considered earlier. Recall for this example: 


0.5 tke 03 0.5 0.2 
A=|0 08 0.2) B=|04 04 02| 2z={0.9, 0.1, 0}, 
(eon al 0.6 03 0.1 
O=fl, 2, 2,3}. 


Thus, using eq (2.16) leads to: 
a(1)=7,b(0, )=0.9-5,(1)=0.9 -0.3=0.27 
a(2)=72,b,(O, )=0.1-b,(1)=0.1-0.4 =0.04 
a,(3)=7,b,(0, )=0-b,(1)=0. 
Next, using the recursion formula (2.17) leads to the following values for 


the forward variable (1) for t=1, 2,..., T-1 and 1Si<N. For example, 


a, (1) = (1a, (7)a,)-2,(0;) 


»& (1)a,,)-B, (2) 


= ( 

= (az, (1)-a,, +Q, (2)-a,, +a, (3)- a5, )-0.5 
= (0.27-0.5+0.04-0+0)-0.5 

= 0.0675 
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a,(2)= (dm (i)a, )+b, (0, Ne 


7 (Ya (1)a,, y): b, (2) 


= (a, (1): A) we Qa, (2). a4 at Q, (3): a3, ) ‘0.4 
= (0.27-0.3+ 0.04-0.8+0)-0.4 
= 0.0452 


a, (3) = (Ya (i)a, )+b; (0, ) 


= (a (i)a;,)-b; (2) 


= (a, (1)-a,, + Q, (2)-ay, + Oy (3)-a5,)-0.3 
= (0.27-0.2 + 0.04-0.2+0)-0.3 
= 0.0186. 


All values of the forward variable are shown below in the matrix Agy: 


0.27 0.04 0 

| 0.06 0.0452 0.0186 
“| 0.016875 0.022564 0.012342 
0.0016875 0.00462274 0.00202298 


which is of dimension: 4x3 (TXN). Similarly, we compute the backward vanable, 
according to eqs (2.20) & (2.21) starting from t=T. Recall from eq (2.20) that: 


Bii)=1, Wi=1,2,3. 


Next, using the recursion formula given in eq (2.21) to compute [3(i) leads to: 
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B; (= Y4,b (0. Bad) 


I 
Me 


2, % ;?, (O, )B (1) 


fee, 
II 
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lI 
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Qj b (3 )8,0)) 


lf 
se 


ry 


11D, (3) + ayy +B, (3) + a3 +B, (3)]- BC) 
=(0.5-0.2+0.3-0.2+0.2-0.1)+1 
alee 


a, 


ll 
© 


All values for the backward variable B,(i) are shown below in the matrix Byy, as: 


0.027582 0.022152 0.009 
0.0726 0.0636 0.03 
0.18 0.18 0.1 
] 1 1 


a 


Note, that the values in the last row are always 1, according to eq (2.20). At this point, 


we can evaluate Pr(O|A), using the forward variable, according to eq (2.18), which leads 
to: 
(OlA)=¥" a,(i)= Yai) =a, (1) + &, (2) + @, (3) + (4) 


= 0.0016875 + 0.00462274 + 0.00202298 
= 0.00833322000000. 


pap Problem 2: Optimal State Estimation 
Problem 2 deals with finding the optimal state sequence I for a given observation 
sequence O. One possible solution to this problem is to maximize the expected number 
of correct individual states by choosing the states i, that are more likely to occur [3]. This 
computation uses the variable y,(i), defined as the probability of being in state qj at time t, 


given the observation O and the model A [3]: 


a aera (2.23) 


yap 


¥:(1) can be expressed as: 


“sem, 


\_ 4,0) 6, (i) 
Ai) = Pr(O | A) > (2.24) 


where ,(t), B,(t) are the forward and backward variables, defined earlier. Replacing eq 


(2.23) into eq (2.14) leads to: 


a,(1)p, (i 
y= ORO (2.25) 
AG) AG. 
i=} 
Note that the normalization factor Pr(O|A) in the denominator of eq (2.24) is needed to 
make y,(i) a conditional probability, which leads to the following constrain being 


satisfied: 
N 
DY A=. (2.26) 
t=] 


Finally, the optimal state sequence i; can be obtained by: 


i =argmax|y,(i)], 1<¢<T. (2.27) 


1<i<N 
A more efficient approach to compute the optimal state sequence uses a decoding based 
on dynamic programming called Viterbi algorithm and shown in Table 2.4. The Viterbi 
algorithm is designed to find the best path (sequence) which maximizes the probability 


Pr(O, IA) [3]. 


ZS 


Table 2.4 Viterbi Algorithm 


Initialization: 


0,() =2,b,(0,), 1SiSN 
9,(i)=0 


Recursion: 


t 


5,(j)= max|6,.,(i)a, b,(0,), for2st<T, 1s j<N 


g,(j)=arg max|6,_, (i)a, 


y 
Termination: 


P* = max|6, ()] 


i> = arg max[6, (i)] 


l<isN 
Path (State Sequence) backtracking: 


i =@,,,(@,,), Fort=T-1,T-2,--41 





Application: 


As an application, let’s go back to example 3 defined earlier in section 3 to find 
the optimal sequence i,, given the model A and the observation O. Recall that, for the 


given model: 


0.5 03 02 O3mOS 02 
A= a5 =| Om Oroee neo 4b), 0 aan OP alee = (ONO) 4122 
0 Cees! 06 03 01 


Thus, according to Table 2.4, we get: 


J (j=EOOn 11a 
6, (1) = 2,b,(O, )=0.9-0.3 =0.27 
6,(2) = 2,b,(O, ) =0.1-0.4 = 0.04 
6,(3)= 2,b,(O, ) =0 
a us 0 
Re (ia, ,( (O,), Pepsin SO: 


a a (ia, b,(o O,) 
7 fe axl 5, (i)a,, b,(O, )=max[6, (1)a,,,4,(2)an,,6, (3)a3, B, 


Similarly, after computing all 6 values, we define the matrix A, such as: 


0.27 0.04 0 

0.0675 0.0324 0.0162 
0.016875 0.010368 0.00486 
0.0016875 0.00165888 0.000486 


A={6,(7)}= 


Then, we compute the @ variable, defined in Table 2.3 as: 


g, (j)= arg max|65,_, (7a, |. 


1sisN 


For example, @2(1) is defined as: 


g,(1)= arg max|6, (i)a,, |= arg max(6, (1 1)a,,,6 2 a 6, (3)as,] 


lsisN 


= arg max(0.27 -0.5,0.04-0,0]=1. 


Similarly, we compute all other ¢,(j) variable, and define the ® matrix: 


© =19,(i)}= 


—=— = — 
No Nw — O&O 
WwW WwW — © 


Next, we find the last state at time T=4 defined as: 


i; = arg max|6, (i)) 


lsisN 


which leads to: 
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i, = arg max (6, (i)| = arg max[6, (1), 5, (2), 6, (3)] 
= arg max(0, (1), 6; (2), 6; (3)]= arg max|0.0016875 0.0016588 0.000486] 
=1. 


Finally, the other state sequence are defined using the backtracking method of the Viterbi 


algorithm given in Table 2.4, which leads to: 


=O (i", ) oma oe 1. 
i; = 9,(i;)=9,(1)=1 
i, =9;(i3)=9,(1)=1 
i; = 9,(i; )=9,(1)=1. 


Therefore, the optimal state sequence for this model and observation sequence is given by 
i={1,1,1,1}, which shows that all balls come from the first urn. Note that repeating the 
same experiment with the observation sequence, O={1,2,3,3}, results in 1={1,2,2,2}. 
Further, note that we always expect a forward moving state sequence, since our model is 
a left-right model. In addition, we can also use the variable 6,(j) to evaluate Pr(OJA), and 
thus we introduce the Viterbi probability Pr,, such as: 


Bie max}6, (i)}. (2.28) 


Using eg (2.28) leads to: 


Pr, = max{6, (i)$= max{6, (1) ,6,(2),6,(3)} 


= max{0.0016875 0.00165888 0.000486} 
= 0.0016875, 


which is consistent (meaning in the same range) to the value 0.00833 found earlier using 
the Baum-Welch probability eq (2.18). Note that we can’t expect to find the same exact 
value (since theoretically is not the same), but we can use both probabilities to reconfirm 


the decision. 
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3: Problem 3: Re-estimation of Model Parameters 


The last problem deals with adjusting the model parameters (A, B, 7), so that 
probability Pr(O|A) is maximum. To that end, we define the variable &(i,j) as the 
probability of a path being in state q; at time t and making a transition to state q; at time 


t+1, given the observation sequence and the model, such as [3]: 


) a a, (i)a,,b, es Boa (j) 


Pr(O [A) ee 


los 


Recalling the definition of y,(i) given earlier in eq (2.25) as the the probability of being in 
state qj at time t, given the observation O and the model A, and using eqs (2.23) and 


(2.24), we can relate y;,(i) to &(i), by summing €(i) over all states j, which leads to: 


y,(i)= SE ey): (2.30) 


j=} 


Similarly, summing y,(i) and &(i) for all t’s, leads to the following result[3]: 


T-1 
yy y, (i) = Expected number of transitions made from state q, , (2.31) 
t=] 


and: 


it 
é,(i,j) = Expected number of transitions made from state q, to state q pe (2.22) 


t=] 
Finally, using eqs (2.31), (2.32) and the definitions of the model parameters, we can 


reestimate the model, according to the Baum-Welch formulas: 


l.2,=y7,(i), 1SiSN, (2.33) 


2 


2.a, = HL 
y Pa ’ (2.34) 
y, (i) 
‘=] 
dv () 
3. bj(k)=-=— (2.35) 
AG) 


We then continue reestimating our model (applying the new model to the variables y and 


€), with the re-estimations formulas defined above, until we reach convergence, i.e., so 
that: 


A iz Acid ° 


new 


Application: Um-ball example. 


Next, we apply the re-estimation formulas to our urn-ball model described earlier: 


Recall: 
0.5 03 0.2 0.3 0.5 0.2 
A= jay;=| 0 08 02198 =(b.;=|04 04 02), 7 =17, f= (09) 0mm 
0 0 1 0.6 03 0.1 


For the observation: O =/1, 2,2,3/, anda single iteration we get: 


0.6256 0.2945 0.0079 0.4362 0.4649 0.0988 
A={aj\=} 0 0.8984 0.1015|, B={b«e}=|0.0712 0.5573 0.37141, 
0 0 1 0 0.4697 0.5302 


II = {u}= {0.8936, 0.1063 0}. 
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Note that, 1) the summation all rows of the matrixes A, B and vector 7% remains equal to 


"wee, 


1, and 2) A is still an upper tnangular matrix, since our model is let-right. Re-estimating 


the model for 10 iterations leads to: 


eel 0 rt 60 0 
A=|0 0.3341 0.6658), B=|0 1 0 | 

0 0 1 0 0.3325 0.6674 
TI ={1, 0, 0}. 


As we can see from the above results, the matrix A still remains an upper 
triangular, which means that the model is still a left-rght model. Examining the B matrix, 
we conclude that, if we are in the first state, we will observe only the 1** symbol (1* row, 


1“ column =1), etc. Finally, the = matrix, shows that procedure starts from the first state. 


F. SCALING 

Numerical implementations of the forward-backward, Baum-Welch or Viterbi 
algorithm may lead to underflow problems due to the small numbers involved in the 
required computations. In addition, the problem may become worse, as the matrix 
dimensions involved in the computations increase. As a result, scaling is required to 
avoid such mathematical problems [4]. The basic idea behind the scaling procedure is to 


multiply the forward and backward variables o,(i) and B,{i) by a coefficient so that the 
scaled G (i) and £,(i) are kept within the dynamic range of the computer. 


Note that we can rewrite the reestimation formula eq (2.34), in terms of the the 


forward and backward variables, as [4]: 


To! 


= ye (i)a,, (0,,, Noes (j) 


ay = (2.36) 


S¥.a,(i)a,b,(O)Bar(J) 


t=] j=l 


Consider the forward algorithm used to compute of the forward variable o,(1), discussed 


earlier. O,(i) is multiplied by the “scaling” factor c, defined as: 


] 


»@ (i) 





Cc, = 


leading to the scaled forward variable defined as: 
G, (i) = ¢,@, (7). 


The scaled coefficient @ (i) can be shown to be equal to [4, pg. 272] : 


N 
Dea) a,b, ( 


=1 


Bane, 


Al 
Le a, -.{ Ja , b,( oF 


i= 


T 


We also define a modified forward variable, as: 
G, (i) = ya &,_,(j)a,b,(O,). 


By induction @,_,(j) can be written in terms of &_,(j) as: 


Thus: 





Sy 4 2 fT os },6,0. j(0,) ae 


1 


The scaled backward variable £ 


i 


(i) is defined as: 
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(2.3%) 


(2.38) 


(239) 


(2.40) 


(2.41) 


(2.42) 


(2.43) 


At this point, we can rewrite the reestimation formula given in eq (2.36), using the scaled 


forward variable such as: 


al 


2,4,{ t a; b,( on ew) 
{= 


Sia AT iN i (2.44) 
2 24| Qa, ai b,( OF, Bb. (i) 
t=] j=l 
Similarly the re-estimation formula for bj(k) eq (2.35) becomes: 
T-1 . pe 
a, (i)B,(4) 
= t=l 
b ;(k) = =. (2.45) 


— 


(= 


Finally, it can be proved that the probability Pr(O|A) and the Viterbi probability Pr, can 


be computed using the scaling factor [4],which leads to: 





Pr(O | A)= e (2.46) 
I<. 
Or, 
log[Pr(O | A)|= -Yiloge (2.47) 
log(Pr, )= max|¢, (i). (2.48) 


G. MULTIPLE OBSERVATION SEQUENCES 


In real word applications we need to train the HMM using multiple trials, to 
ensure robustness in the recognition/classification process. Each training trial signal 


produces an observation sequence. Assume that there are K trials, i.e., K observation 


sequences. We denote the set of observation sequences O such as: 
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oF io 0), wu 0}. (2.49) 
where: 

9 &) = 10),0®),.-,0, (2.50) 
is the k” observation sequence. Provided that each observation sequence is independent 
of each other and identically distributed, we want to find the model which maximizes the 


probability: 


Pr(O | A) = [[Pr(o® |A)= ee | (2.51) 


k=1 k=1 


Therefore, the multiple observation reestimation formulas using the scaled variables are: 


S_. Aha (kK) \ Qk (; 
ies De (i a,b, (0! t+] (j) 
k 


Fp Sec (2.52) 


b, = = a ee (2,53) 


Note that we don’t reestimate 7, so we keep it unchanged through the iterations. 


a2 


II. ISOLATED WORD RECOGNITION 


This section presents the application of HMMs to speech recognition. 
Specifically, we show how the model parameters have to be selected and adjusted using a 
limited 4-word recognition example. Our goal is to show how a more genenic classifier 
can be set-up through the speech recognition example. In addition, we point out the 
potential difficulties one may encounter while setting up and implementing the software 


for such an application, due to computer numerical precision limitations. 


A. GENERAL HMM TRAINING AND TEST PROCEDURE 


This section describes the application of HMMs to classification in the context of 
speech recognition. The overall procedure ts illustrated in Figure 3.1. First, labeled data is 
used to train the model. These specific training signals may be multiple trials of the same 
type of signal, 1.e., belong to the same class, or belong to different classes. For example, 
multiple trials of the same word may belong to the same signal class in speech 
recognition applications. In such a case, all words used in the recognition set-up 
constitute the dictionary. Next, information uniquely characterizing each class needs to 
be extracted from the signal classes. Thus, each signal is split into T segments, and some 
useful features extracted from each. Feature vectors may include LCP or cepstral 
coefficients, energy, etc... Thus, the initial signals are converted into a set of continued- 
valued vectors. Next, this set gets converted into a sequence of discrete vectors using 
vector quantization (VQ) [4,5], which will be described further in Section 2. The set of 
M discrete symbols forms the codebook. For example, assume we have a speech signal 
divided into 4 segments (i.e., T=4), and that the set of symbols representing one signal 


segment is a number ranging from | to 8. Thus, a possible observation sequence of the 
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k™ signal may be given as O* = {7, 3, 4, 8}. At this point, we can check whether the 
features extraction method and dimension of the codebook M makes sense by comparing 
the observation sequences obtained for the signals of the same class, as it is reasonable to 
expect some similarity between the resulting sequences. 

Recall that the set of observations derived from each class of training signals is 
the only information used to train a class-specific model A{A, B, 7}. First, an initial 
estimate of the model is required to apply the Baum-Welch algorithm, as described 
earlier in Section LU. Initial values may be set randomly so that they satisfy the model 
constraints and converge to the correct type of model, 1.e., the initial matrix A 1s to be 
selected as upper triangular for left-to-nght models so that it converges to an upper 
triangular matnix. The Baum-Welch algorithm iteratively estimates the model and stops 


when there is no significant model parameter changes between two successive iterations. 
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Figure 3.1 HMM General Training and Testing Procedure 


a) 


Unlabeled signals are used during the testing phase to identify which class they belong to. 
First, testing data are pre-processed to generate the codebook vectors and corresponding 
observation sequences, following the same process as that used during the training phase. 
Next, the set of the probabilities observing the tested signal O, given the Kk" model, 
Pr(O|A™) is computed using either the Baum-Welsh or the Viterbi algorithm for all k™ 
models, as discussed earlier in Section II. Finally, the model type is selected by choosing 
that with the highest probability, Pr(O|A") . 

Next, we describe a simple 4-word recognizer designed to recognize the words: 
Statistics, Microsoft, Instructor, and Professor. Three trial words are used for training and 
one word is use for testing for each class. 


1. Data Creation and Preparation 


All data were recorded on the same machine running Windows-98 Sound 
Recorder with a sampling rate of 8000Hz sampling and 8bit mono encoding. One male 
speaker was used. However, note that there are some variations between the trials so that 
no word is pronounced twice exactly the same way. An energy detector was applied to 
remove silence before and after each word, resulting in word of about 9000 points. Next, 
each signal was interpolated to 10000 points to obtain trials with the same length. Finally, 
each data was sent through a pre-emphasis filter with transfer function H(z)=1-az" with 
a=0.98, to emphasize the relative energy of the high-frequency spectrum which contain 
useful information. The MATLAB software implementation is presented Appendix A.3. 


Ze Feature Extraction 


We varied the length and number of segments (time frame), and the specific type of 


feature parameters until we obtained a combination which lead to correct classification. 
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The final set of parameters was derived by each word into T=7 segments using a 
rectangular window, with an overlap a 10%, where the time length of each trial was 
about lsec. The following eight parameters were extracted for later use in VQ from each 
segment: LPC coefficients from a 7" order filter derived with the covariance method [1], 
and the energy of the section. Finally, we normalize the value of the energy dividing all 


values by the max energy. 


Sh Vector Quantization 


Vector quantization (VQ) is a scheme, which maps a sequence of continuous- 
valued vectors into a sequence with a given number of discrete vectors, called the 
codebook [7]. Therefore VQ can be viewed as some type of encoding scheme, where the 
encoder Y assigns a channel symbol y(x) from an ensemble of M symbols to each input 
vector x={Xo, X1,..-, Xk1} [7]. Note that there is no need to define a decoder, as the 
discrete sets of parameters never get translated back into the original vector. 

Basically, VQ partitions the set of coefficients into M disjoint sets. Each set is 
represented by a single vector {vm}1l<m<M, which is a centroid of the vectors in the 
coefficient set assigned to the m™ region [4]. 

Note that there is a distortion penalty associated with VQ, as all feature vectors 
are represented using a set of M codebook vectors. The larger the dimension of the 
codebook, the smaller is the overall distortion between orginal feature vectors x and 
codebook vectors x;. The distortion measure d(x,x,) resulting from the codebook selection 


process can be represented using the Euclidean norm as [7]: 


(3.4) 








een |x — x. : 


aii! 


For our purpose two schemes were applied to derive the codebook vectors: 1) a 
competitive Neural Networks implementation, and 2) a K-means (or LGB algorithm) 


algorithm [7, 8]. 


a) Competitive Neural Network Implementation 


Input Competitive Layer 





Figure 3.2 Competitive Neural Network 


A competitive unsupervised neural network (NN) implementation was 
selected to compute the codebook, as shown in Figure 3.2 [2]. Basically, this NN can be 
viewed as a Clustering scheme, where weights associated to each neuron are used to 
assign a symbol M to each word segment. The software implementation is presented in 


Appendix A.1. 
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Figure 3.3 Diagram of coefficient vectors and VQ vectors, Euclidean distance for the 
word “statistics”, using a competitive neural network 


Figure 3.3 shows seven original feature vectors, and the corresponding 
quantized vectors obtained with the competitive NN. Training took about 2000 epochs 
and 50mn with a Pentium-III 450MHz. Each initial feature vector is mapped into one of 
the M=4 quantized vectors. For example, vectors 2 and 3 get mapped to the same fourth 
class due to their consistency. Similarly, vectors 5, 6, and 7 get mapped to the 3™ class of 


the quantized vectors. 


b) K-means Scheme 


The second method considered is the K-means algorithm, also called the 


Lloyd or Linde-Buzo-Gray (LBG) algorithm [9]. The software implementation is given 
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in Appendices A.13 and A.14 [8]. Basically, the K-means algorithm is a clustering 
scheme, which iteratively finds a set of k, quantized vectors into which all training 
vectors get mapped to with minimum distortion. The number of clusters increases 
iteratively by splitting the existing quantized vectors obtained at each iteration, until a 
desired number of quantized vectors or distortion levels is obtained [9]. Thus, the size of 

the codebook is a power of two, 1.e., 2, 4, 8, 16, etc... The LGB algorithm is illustrated 
for the word “Microsoft” in Figures 3.3 and 3.4 for codebook sizes equal to 4 and 32 


respectively. 
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Figure 3.4 Diagram of coefficient vectors and VQ vectors, Euclidean Distance for word 
“Microsoft,” using the K-means (LGB) algorithm and M=372. 


Note the distortions are much smaller when M=32 than when M=4. This 
is to be expected as a larger codebook size allows more flexibility. However, a large 
codebook may not be always desirable as it may lead problems in evaluating and 
comparing the resulting probabilities Pr(OlA), as we will see later. The K-means 
implementation is a little faster than the NN technique, however it is restricted to 
codebook sizes which are powers of two. No such restriction is needed for the NN 
implementation. 

Note that an initial selection of the number of segments T and the size of 


the codebook M may be made by observing a small selection of the resulting quantized 
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vectors. The basic idea is to select a combination of parameters, which lead to some type 
of consistency in the resulting quantized vectors in a given class. 


4. | HMMTraining 


This section describes the application of the HMM training theory presented in 
Section II to the speech recognition example. The HMM implementation uses MATLAB 
5.3, and is presented in Appendix A. MATLAB may not be the most desirable software 
language as it is relatively slow, however it was selected for ease of implementation to 


test the concepts. 


a) HMM Initial Conditions 


Usually, a left-to-right model is preferred over the more general ergodic 
one in speech applications because model states can be associated with time in a 
straightforward manner [4]. Thus, we selected a left-to-right model, where the matrix A 
is upper triangular. In addition, we also prefer the model to start from an early state so 
that it can pass from all possible states. As a result, we didn’t use random values for the 
vector 7, but instead we initialized the first coordinate to be quite large, the second one 
smaller, and so on, and applied the constrain that the summation of all values in 7 1s 
equal to one. We use random values satisfying the appropriate constraint for the matrix 


B. 


b) Number of States N 


HMM states are called “hidden” because all information about the state 
sequence is not accessible directly, since the only information available is given by the 
observations. According to Rabiner [4], there are two schools of thought for the physical 


meaning of the number of states; the first one states that it represents the number of 
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sounds (i.e., phonemes) for each word, which is usually a number between 2 and 10. The 
second interpretation is that it represents the average number of observations in a spoken 
version of the word. Basically, we can assume that the number of states represents the 
number of distinct sounds (e.g., phonemes or syllables) of the word. In our case a 
number of states N equal to 4 was deemed appropriate, since our vocabulary contains 
only 4 words. 

Further investigations may be required when applying HMMs to other 
types of signals if we cannot relate the number of states to some physical meaning behind 
the data. At worse, the number of states can be selected by trial and error at that leading 
to the highest recognition results. Note that selecting too small a number of states may 
prevent from differentiating between classes. Simulations showed that selecting too large 
a number of states may result in overly long training time and numerical instability in the 


implementation which cannot be controlled with scaling. 


Cc) HMM Re-estimation 


The HMM re-estimation iterative scheme implemented uses the multiple 
observation sequence technique described earlier in Section 2. The MATLAB software 
implementation is given in Appendix A.7 to A.9. First, we re-estimate the model 
parameters for every different trial, and then we apply these results in the re-estimation 
formula in eqs (2.52) and (2.53). 

The HMM re-estimation step uses the scaled forward and backward 
variables to avoid numerical instabilities, and its implementation is given in Appendix 
A.6. In addition, we ran all MATLAB files using a scaled fixed point format with 15 


digits (format long). Since, some computations still resulted in “divided by zero” errors 
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after applying these corrective measures. We forced the denominator quantities to be 
equal to the smallest floating point number whenever there were found to be smaller. 

We used multiple iterations for each single observation HMM re- 
estimation step until convergence was reached. We noticed that the model doesn’t 
necessarily converge to the same parameters for the same observation, when repeating 
the re-estimation procedure with random initial conditions. However, correct decision is 


still achieved. 


d) Scoring 

The system is ready to test any word and categorize it as one the four class 
types after the four models rm k={1,2,3,4} derived separately for each word (Microsoft, 
Statistics, Instructor and Professor) are identified. We repeat the same feature extraction 
scheme for the testing words, using either the neural network or the codebook derived in 
the K-means algorithm during the training process. Let’s assume that the observation of 
the testing word is O. During testing, we evaluated the probability Pry(OlA™) for 
k=1,...4, using either the Baum-Welch or the Viterbi algorithm, and selected the model 


which gives the highest probability Pr(OlA™) (Appendix A.10, A.11). 


Classification of testing words using all models with parameters: 
M=4, N=8, T=7 
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Tested words: 1: Microsoft, 2: Statistics, 3: Instructor, 4: Professor 
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<i : Correct Decision (100% success) 


Figure 3.5 Scoring of all 4 testing words (1: Microsoft, 2: Statistics, 3: Instructor, 4: 
Professor) given each model X) Parameters: M=4, N=2, T=7. This system 
performed 100% successful decisions. 

Figure 3.5 presents the results obtained when the number of segments T is 
equal to 4, the number of symbols M is equal to 4 ( computed with the LGB algorithm), 
and the number of states N is equal to 8. The four horizontal bars contained in each word 
represent the probabilities 10logio(P(O|A)), k=1,...4. Thus, the highest probability is 
that closer to the OdB point, which represents the point of probability 1. The left column 


plots represent the results obtained with the Baum-Welsh algorithm while the nght colum 


plots represent those obtained with the Viterbi algorithm. Recall that the four models 


45 


MP 72) 7 and A were trained with three trials for every class. The software 


implementation for this testing step is given in Appendix A.11. Results show that correct 


classification is obtained in all 32 cases, as shown in Figure 3.5. 
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Figure 3.6 Scoring obtained for all four testing words (1: Microsoft, 2: Statistics, 3: 
Instructor, 4: Professor) given each model x. Model Parameters selected: M=4, 
N=2, T=7. This HMM has too few states (N=2) and cannot classify the word 
‘instructor’ correctly. 
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Figure 3.6 presents the results obtained when selecting too small a number 


—, 


of states (N=2). All other parameters are kept the same as those in Figure 3.6. Note that 
the algorithm reaches the wrong decision for the word “instructor” which was found to be 


“professor.”’ 


B. CONCLUSIONS 


This section presented an HMM-based classifier applied to a simple 4-word 
recognition problem. We implemented and tested the classifier using MATLAB, and 
described how we selected the various parameters which need to be selected to set-up the 
classifier. Next Chapter considers the application of the HMM-based classifier to 
seismo-acoustic signals to differentiate between two types of mine-like objects buried in 


shallow sand. 
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IV. _HMM-BASED CLASSIFICATION OF SEISMO-ACOUSTIC MINE 
SIGNALS 


HMMs can be applied to various types of classification problems. This Section 
presents the results obtained for the classification of mine-like objects buried in sand. The 
mine data was obtained from Prof. Muir who heads a buried mine detection project 
started in November 1996 at the Naval Postgraduate School. Initial results obtained are 
described in Gagham [10], Fitzpatrick [11], and Hall [12]. 

The NPS mine project is a continuation of work started at the Applied Research 
Laboratory of the University of Texas at Austin, and is sponsored by the Office of Naval 
Research. The main goal of the project is to study the development of a seismo-acoustic 
sonar for the detection of buried ordnance using guided, seismic interface waves. Earlier 
ARL and NPS results showed that seismic interface (Rayleigh) waves can be used to 
detect mine-like objects buried in sand [10-12]. The goal of the current NPS program is 
to develop an improved seismic source to evaluate the feasibility of using a seismo- 
acoustic sonar to detect buned ordnance in the beach and surf zones. The seismic waves 
were generated by two actuators, and were measured by two three-axis sensors- 
geophones. Buried, mine-like objects, ranging from 71kg to 290kg, and at ranges of up 
to 5 meters were echo-located by applying a basic polarization filtering signal processing 
scheme. 

This section applies the HMM concepts derived earlier to distinguish between two 
types of mine-like objects. Two types of experiments were conducted: 1) classification of 


two different mine-like objects at the same range with multiple weights for each mine 
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type; and 2) classification of two different mine-lhke objects at 3 different ranges with the 


same weight for each mine type. 


A. BASIC EXPERIMENT INFORMATION 


This section does not present any theoretical, physical or experimental details for 
this project, as they may be found in in [10-12]. However, we do provide some basic 
information regarding the physical nature of the signals under study. 

The beach site used for collecting the data is a stretch of U.S. Navy-owned beach 
directly seaward of NPS, Monterey, CA, and shown in Appendix B.2. This area 
measured roughly 150 feet in length running parallel] to the waterline and vaned from 20 
to 50 feet from the high-to-low water mark. Generally, the sand conditions varied due to 
the different waterline distance, resulting in changes on some of the sand characteristics, 
such as density and moisture in multiple depth layers. 

The equipment configuration for the mine detection is illustrated in Appendix 
B.6. Basically, two actuators produce the Rayleigh waves described in [11]. Rayleigh 
waves have distinctive features that make them identifiable in a complex seismo-acoustic 
wavefield. The most important property of Rayleigh waves is the elliptical particle 
motion produced by their passage. Their unique characteristic is the 90° phase shift 
between their horizontal] and vertical components, which results in elliptical motion, as 
illustrated in Appendix B.7. Thus, the placement and operation of the actuators, as 
illustrated in Figure 4.1, can be categorized into two different configurations: 1) actuators 
operated horizontally; and 2) actuators operated vertically, with a 90° relative phase 


difference between their drive signals [11]. 
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Figure 4.1 Actuator placement for Rayleigh waves generation [11]. 


The generated seismo-acoustic wave travels through the geophones and meets the 
target, as illustrated in Appendices B.6 and B.8. Next, the reflected wave passes through 
the geophones again, and the received signal is used to detect the target presence. The 
signal is collected, as described in [12], and the signals for the relative x, y, and z-motion 
are shown in Appendix B.5. 

Note that the target reflection is not visible from the raw signals. Thus, vector 
polarization (VP) [12] was used to extract the Rayleigh waves from the raw information. 
The VP step takes advantage of the 90° phase shift between vertical and radial 
components to pull the target information out. This is accomplished by applying the 
Hilbert Transform to the radial-vertical signal pair [11], the complex crossed-power is 


computed such as [12]: 


Po =Nhitten *V hitbert (4.1) 


the imaginary component of P,, is essentially proportional to the intensity of the seismic 
wave due to the 90° phase shift between vertical and radial components [12]. The 
polarity of the imaginary power is associated with the rotation of the elliptical motion of 


the wave (e.g., a negative value corresponds to an anticlock-wise motion). Real and 


Si 


imaginary components of the cross-power Py, are illustrated in Appendix B.9. Appendix 
B.8 shows the incident wave passing the geophone at 10ft. The reflected target signal 
Strength is quite small, and we need to expand the scale to see it, as explained further 
later. 

Two mine-like objects were used: 1) a cylinder weighing 150]b (68kg), 5ft long, 
8in (20cm) in diameter, with a %4-in (0.6cm) wall thickness; and 2) a U.S. Navy power 
can or “power keg” with the shape of a cylindrical sheet metal can 18in (46cm) high and 
24in (61cm) in diameter, weighing 16lb (7kg). These objects are shown in Appendix 
B.3. Additional weights were used to vary the objects total weight. These additional 
weights made the maximum weight of the cylinder and the keg equal to 618lb (280kg), 
and 640lb (290kg) respectively. The cylinder was always buried on its side with the 
cylindrical axis horizontal and in the direction of the wave propagation. The powder keg 
was always buried upright, with the cylindrical axis vertical to the direction of the wave 


propagation. Further details are given in Appendix B.4. 


B. SIGNAL SELECTION 


Recall that the target signal is visible after processing the raw information using 
VP. Thus, we use the signals obtained from the imaginary portion of the cross power for 
our application only. We also focus on the experiments conducted on the 6" and 10" of 
November 1998, were the cylinder and the keg targets were used, as shown in Appendix 
B.6. Two parameters were varied during these two days in addition to the sea conditions: 
target weights and distances from the geophones. Ten tnals were conducted for each 
experiment. The distances between the different pieces of the equipment (actuators, 


geophones, target) were set as shown in Appendix B.6. As a result, the actuator- 


geophone distance was set at 10ft, the actuator-target distance: 16ft, and the actuator- 
target-geophone (reflected) was 22ft Fa both targets (cylinder and keg). Appendix B.8 
shows the average imaginary cross power obtained for all trials associated with a specific 
target and weight, scaled near the distance of 22ft. Appendix B.9 shows that the reflected 
target signal becomes stronger as the target mass increases. Appendices B.10 and B.11 
plot the imaginary cross-power signals obtained from the cylinder and the keg targets 
respectively. We used the signal from the 2"° geophone for the cylinder case, as it is the 
strongest of the two, and the signal from the 1 geophone for the powder keg, as the Des 
geophone doesn’t show any received signal. However, the signal measured at the oe 


geophone for the keg experiment is used as an example of background noise (no target) 


signal. This background noise is used as part of the testing data, as described later. 


Cc SIGNAL PREPARATION 


Plots contained in Appendix B.10 and B.11 show there is some consistency 
between the signals obtained from multiple weights of the same target. Thus, we elected 
to consider all signals associated with a given target to the same class, independent of the 
weight. Five trials were used for training, and the 6" trial was used for testing. The 
following weights were used for the keg target: 1) 224lb; 2) 432lb; and 3) 536lb. 
Similarly, the following weights were used for the cylinder target: 1) 3641b; 2) 468lb 3) 
572lb. Next, we assumed that we detected the presence of a non-background signal 
using some type of threshold detector, and truncated the signals used for the HMM set-up 
around the position of the target location. Signals used for the HMMs training and 
classification are shown in Figure 4.2. Note that it would be useless to include the 


incident signal received at a geophone from the actuators, as it doesn’t contain any 


ae 


information about the target, and is much stronger than the signal reflected from the 
target. Each signal is 400 data points long. Two classes (keg and cylinder) are 
generated. At this point, the goal 1s to recognize the shape of the target, and each class 


needs to be characterized by a set of features, as explained next. 


Powder Keg 224!b Cylinder 364|b 
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Figure 4.2 Imaginary component cross-power for the keg and cylinder targets for 3 
different weights, truncated near the target location. Each plot illustrates 6 different trials 
(for the same target, and weight) 


D. FEATURES EXTRACTION/VECTOR QUANTIZATION 


Note in Figure 4.2 that the signals contain little spectral information. The 
frequency of the Rayleigh waves obtained during the November 6" and 10” experiments 
was 80Hz. As a result, the signals contain very little useful spectral information. We 
tried to apply LPC feature extraction, and other similar spectral coefficients, but we 
found little or no consistency between the trials of a given class. Therefore, we simply 


used the time-domain information directly. We segmented the signals into two segments 
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(T=2), and decimated each frame by a factor of 20:1, resulting in 10 data points over each 
segment, as shown in Appendix C.4. erally. we normalized those points to compensate 
for the differences in signal strength coming from different weights, by dividing every 
value with the maximum one. 


We applied the LGB algorithm at the VQ stage and selected M=8 symbols. The 


code implementation is given in Appendix C.3. 


E. HMM TRAINING FOR THE MULTIPLE WEIGHTS EXPERIMENT 


HMMs need data to be trained. However, in practice the availability of data may 
be seriously limited for various reasons and cross-validation needs to be used [13]. 

This study has only a limited number of trials available: six trials for every type of 
mine-like object. Thus, we used cross-validation by successively selecting five trials for 
training and the last one for testing. The procedure was repeated six times rotating the 
testing signal each time. The code implementation is given in Appendices C.5 and C.6. 
We found that the best performance is achieved with the number of HMM states N equal 
to four. 

We created two ergodic HMM models: one for the keg class and one for the 
cylinder. Thus, we used a total number of 5(trials)x3(types of weight)=15 signals for the 
HMM training for every class (i.e., model). We selected an ergodic model as it 


performed better than the left-to-nght model for the data. 


F. MULTIPLE WEIGHTS SCORING AND RESULTS 


The MATLAB file sequence.m given in Appendix C.6 performs all HMM 


training and scoring. Note that the testing signal is rotated each time, and the results 


5° 


plotted, as shown in Figure 4.3a (cylinder training case) and in Figure 4.3b (keg training 


case). 
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Figure 4.3a HMM testing results for the cylinder model (multiple weights). 
Prpw(O|Acytinder) and Pry(O|Acyiinder) for all testing signals and all rotations (every row). The 
model Acylinder 1S Created by training all cylinder signals (5x3=15). The decision is correct 
whenever the tested signals (4,5,6) have a higher Probability (.e., closer to 0 on a dB 


scale) than 1, 2, 3, or 7 (keg signals and background). Each row represents one of the six 
iterations of the rotation between testing and training signals. 
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Figure 4.3b HMM testing results for the keg model (multiple weights). 
Prpw(OlAkeg) and Pry(O|Axeg) for all testing signals and all rotations (every row). The 
model Axeg is created by training all keg signals (5x3=15). The decision is correct 
whenever the tested signals (1, 2, 3) have a higher probability (i.e., closer to 0 on a 
dB scale) than 4, 5, 6, or 7 (cylinder signals and background). Each row represents 
one of the six iterations of the rotation between testing and training signals. 
We used as background (non target) signals, the six trials obtained from the 2"° 
geophone during the keg experiment which didn’t track any target signal. Figures 4.3a 
and 4.3b show that an decision error for one case only. The signal considered in that set- 


up is the 4" tested signal, i.e., the cylinder with weight of 368lbs. Figure 4.2 shows that it 


is the weakest signal of all the mine signals. The total number of testing data is 


DT 


2(mines)x3(weights)x6(rotations)=36, thus the overall classification performance of the 


system is (36-1)/36 =97.2%. 


G. MULTIPLE DISTANCES EXPERIMENT 


Up to this point we showed that we can recognize the shape of the mine-like 
objects for a fixed distance. Next, we apply the HMM-based classifier to recognize these 
objects located at three different distances from the actuator, as we need to be able to 
detect and recognize a mine independent of the distance from the sonar equipment for a 
more realistic set-up. Experiments conducted on the 6" and the 10" of November 1998 
provide the following data by moving the location of the geophones between the actuator 
and the target which stays fixed at 16ft, as illustrated in Appendix B.6. Thus, the 


following three set-ups are available: 


e Distance between actuator and geophone equal to 6ft, resulting in the total 
distance (actuator-geophone-target-geophone) equal to 26ft, 


e Distance between actuator and geophone equal to 8ft, resulting in the total 
distance (actuator-geophone-target-geophone) equal to 24ft. 


e Distance between actuator and geophone equal to 10ft, resulting in the total 
distance (actuator-geophone-target-geophone) equal to 22ft. 


Note that the reflected signals obtained for the cylinder for total distances equal to 22ft 
and 24ft were very weak (in fact, same range as that of the background noise). Thus, we 
used this data for the powder keg case only. 

We created two classes again. The first class contains the keg data with weight 
224lbs obtained for the 6ft, 8ft, and LOft experimental set-ups, and 6 trials at each 
distance. The second class is composed of the cylinder data with weight 364lbs for the 
experimental set-up of 10ft. Six trials are also available for the experimental set-up. All 


signals were segmented again around the target location, as done before, resulting in 
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signals with length equal to 400 points, as shown in Figure 4.4. Again, Figure 4.4 shows 
that there is also a viewable Soe for the mine-like keg object. No conclusion can 
be drawn for the cylinder, as the signal was too weak to be usable. 

We used the same feature extraction as that considered earlier for the multiple 
weight experiment (2 segments, decimation, VQ, number of symbols M=8), and the 
implementation 1s presented in Appendix C.9. The model selected is ergodic with four 
States (N=4). Cross-validation was used again due to the limited amount of data 
available. Scoring results obtained for the keg class model and the cylinder class model 
are presented in Figures 4.5a and 4.5b respectively. The MATLAB files implementations 


for the HMM training/testing using multiple distances case are presented in Appendices 


C.10, C.11, and C.12. 





Data Points 


Figure 4.4 Imaginary component of the cross-power for 3 different distances for the keg 
target and one distance for the cylinder. Each plot illustrates 6 different trials (for the 
same target, and distance). 
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Figure 4.5a Prpy(OlAreg) and Pr/(O|Axeg) for all testing signals and all rotations 
(every row). The model Axe, is created by training all keg signals (Sx3=15). The decision 
is correct whenever the tested signals (1, 2, 3) have a higher probability (i.e., closer to 0 
on a dB scale) than 4, or 5 (cylinder signals and background). Each row represents one of 
the six iterations of the rotation between testing and training signals. 
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Figure 4.5b  Prpw(OlAcyiinder) and Pr (OlAcyiinger) for all testing signals and all 
rotations (every row). The model Acyiinder is created by training all cylinder signals 
(5x1=5). The decision is correct whenever the tested signals (4) have a higher probability 
(i.e., closer to O on a dB scale) than 1, 2, 3, or 5 (keg signals and background). Each row 
represents one of the six iterations of the rotation between testing and training signals. 


Figures 4.5a and 4.5b show that the system performs 100% correct detection, 
which seems to indicates that the HMMs took advantage of the consistency of the signals 


at different distances. However, additional] data is needed to make further conclusions. 


H. CONCLUSIONS 
This section described a HMM-based mine-like object classification system using 
the seismo-acoustic waves provided by the NPS project [10-12]. Initial results indicate 


that the HMM-based classifier can recognize the type of mine-like object, independent of 
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the object weight with a 97% accuracy. Results also indicate that it can recognize the 
object type at different distances with a 100% accuracy. However, the experiments were 
conducted with very few data, and further work needs to be done to confirm these initial 


findings by using a larger data set. 
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Ne MINE-LIKE OBJECT RECOGNITION USING NEURAL NETWORKS 


This chapter considers the application of a back-propagation neural network 
classifier to the same data as that considered earlier with the HMM-based classifier. The 
goal of this chapter is to compare the resulting performances obtained with the two 


different implementations. 


A. NEURAL NETWORK DESCRIPTION 


NN input feature vectors are those obtained after the VQ step described earlier in 
Chapter If]. Thus, we have two signal classes again: the keg and the cylinder class. We 
select a back-propagation feed-forward NN (BPNN) with one hidden layer and two 
outputs (one for each class of mine-like objects). Note we do not specific a target output 
for the background signal as: 1) there is no consistency between the different background 
signals available, and 2) we do not use this classification technique to confirm a target 
existence. 

The network specific structure is described in Figure 5.1. BPNNs use the steepest 
descent algorithm, or variants of it, to find the weights which minimize the squared error 
between target and network outputs [2,14]. BPNN may use one or more hidden layers of 


sigmoid neurons, and one output layer of linear neurons[14]. 
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Figure 5.1 Hidden and output layer of neurons of a backpropagation feedforward 
neural network. 


B. MULTIPLE WEIGHTS SET-UP 


We used the observations from the 3 different weights for each mine-like object 
(keg and cylinder), 5 trial each, as input vectors p to train the neural network. We used 
the numbers “1” and “2” as targets for the neural network, such as “1” for the keg class, 
and “2” for the cylinder class (see Appendix D.1 for the MATLAB code). Finally, we 
tested the 6” trial of each signal by computing its output value to the trained NN. The 
test signal is recognized as a “keg” signal when the output is close to 1, and it is 
recognized at a “cylinder” signal when it is close to 2. Note that the NN output is not 
binary, so we set a range around output values 1 and 2, around which the data is said to 


belong to one class or the other. The specific range is set at 3% above or below the target 
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output values 1 and 2. The data is considered as not belonging to either class (other) 


when the NN output value is outside that range. 


We also tested the background signal defined for the HMM-based classifier. We 
rotated training and testing observations six times, as done with the HMM-based 
classifier. Results for the multiple-weights case are given in Table 5.1. Note that the 
system recognizes all testing signals, except the first background signal, which is detected 
as a cylinder-class signal because the associated NN was equal to 1.999 (within the +3% 
range of the target output value 2 associated to the cylinder-class). Thus, the overall 
recognition performance is 97%, as the total number of testing signals is 42 (7signals x 6 


rotations). 
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Table 5.1 Multiple weights NN outputs/decisions © 


Testing Signal 


ffs 0.9999 | 0.9999 | 0.9999 | 1.9999 | 1.9999 1.9999 | 1.9999 


y V V V V V V 


1.0117 geoos Eee. 1.9977 2.0053 2.0053 0.2342 
V: Correct Decision 
X: Wrong Decision 
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C. | MULTIPLE DISTANCES SET-UP 


Similarly, the NN is used for the multiple distance set-up, as considered earlier in 


Section IV. Results are presented in Table 5.2, and the code is given in Appendix D.2. 


Table 5.2 Multiple distances NN outputs/decisions 
Testing Signal 
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Table 5.2 shows that all decisions are correct, except in two cases: 
1) the first trial of the keg data at 8ft gets detected as background, 
2) the first trial of the background signal gets detected as cylinder-class data. 


Thus, the overall classifier performance is 93% for this set-up. 


D. CONCLUSIONS/COMPARISON WITH THE HMM-BASED CLASSIFIER 


These results show that the two classifiers performed in a similar manner. The 
overall performance was 97% for the multiple weights set-up, while the HMM-based 
classifier performance was 100%, and the NN 93% for the multiple distance set-up. Note 
that no background signal was used during the NN training, and that we assigned as 
background data which NN output didn’t fall into the prescribed range. 

Finally, we measure the speed of the execution of the MATLAB file for the 
evaluation of all the results for the multiple distances case. Thus, in a Pentium II 
450MHz, 128MB Ram, the execution time for the HMMs was 5s, and for the NN was 


15s, thus the HMM was 3 times faster that the NN (training and testing). 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


This study presented an introduction to HMMs and their applications to 
classification problems. HMMs require the selection of consistent characteristics between 
the classes to perform well. In addition, many other parameters have to be carefully 
selected to set-up the HMM to have it perform successfully. We implemented the HMM 
using MATLAB and tested it on a simple four isolated word recognition problem. 

The HMM.-based classifier implemented in this study performed very well on the 
limited (two classes) mine-like object recognition experiment. Results also show that its 
performance is similar to that obtained with a BPNN. However, the classifier needs to be 


run using larger data sets to confirm the present findings. 


B. RECOMMENDATIONS 


The reflected signals from the targets were quite weak, especially when the target 
distance was large, and the mine weight low. Furthermore, the received reflected signals 
from the two geophones sometimes had different strengths, (see Appendix B.10), or 
weren’t received by both geophones, (see Appendix B.11). Thus, additional data is 
required to further the classification/recognition process, and data collection is in 
progress, as of September 1999 [15]. Finally, a more complicated HMM-based set-up 
would probably be needed, if we want to train the classifier with different types of mines, 


possibly using a higher number of segments (T), and/or symbols (M), and/or states (N). 
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APPENDIX A. HIDDEN MARKOV MODEL MATLAB PROGRAMS 


This appendix contains the various MATLAB programs used for the HMM-based 


classification of isolated words. 


dA 


APPENDIX Al. FEATURES EXTRACTION/VECTOR QUANTIZATION USING 
NEURAL NETWORKS 


% Filename: training.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Features Analysis (LPC & Energy) of 4 words (Statistics, 
Mie rOSsoOt Co LnStrucEor, 


% and Professor), 3 trials each + 1 test word, vector 
quantization using a competitive layer neural network 

% for HMM classification use. 

clear 

%$ M: dimension of codebook, N: number of states 

M=4;-N=8; 


rand(’seed’,1); 

% the mat files contain the word signals 

$load stat.mat;load micro.mat;load instrprof.mat; 

load words.mat; 

% signals interpolation 

% coeff is the fuction that gives back the LCP + energy coefficients. 
%$ T 1s the number of segments 

[estas Ti=cocftt (stall); |cstaZ, 1] =cocrt (staz)e lcstas, T)=coeft (Stas); [esea 
4,T])=coeff(sta4); 

[emuecnolm hi =CGOetr (m1eCrod)> (eEmiLcroZ, 1] =COcEL(muenoZ)- 

femperes, 1) =coecti (micros) = [emiero4, 1) =ecoctaimiuerod) ; 

feinstnl «| —coe-t (unstrl); leinser2, 1 )—coctrUimstrz) - 

fermstrsy . | =coctm(imstrs) -[ermstrtest, 1] =ecockt (instrtest ) > 

[eproriy Ei —coctt(prorl); (eprerz, T|=CcCoeckr(pror2 

[epror3 ,Ti=Coctt (prot3) >; [eproftest, T|=coett(prottest) ; 


mea 10:0; 

% we convert the results by a constant, because the NN works better 
convert=2000; 
cstal=cstal*convert;cstaZ=csta2*convert;csta3=csta3*convert;cstad=csta4 
FCOnVeEE et; 

Mieco /=—OmEGcrol  COnvenre, cml croZ-Ccmtensea.- GONVere - emi cCro3-Cmicros “Comverse 
emus hoo—-Cm1cro4* Convert: 

Cinserel—-cinstad “Convert, CINSEY2—-cinsers CONVEert ;-CIinStr3=cinst.r3s * converse 
,einstrrest—-Cinstrrest “Convert; 

CDEOLh=Cprot ly Convert, CprOrZ=Ccpreolr2 “Convert ; Cproft3=Cprors “Convert - Coron 
test=cproftest*convert; 


% vector quantization with NN 

pl=(Cstal>cstaz, cstac. 

CMilCrol -Cmicror, Cmlcnos-GCimnst or Cinstry-Cinstrs)]- 
Pia b=Masamasc( pi). 


net=newc (pri,M,1); 
net.trainParam.epochs=1000; 
net.trainParam.show=100; 
net=train(net,pl); 
w=net.IW{1}; 


ys=sim(net,cstal’); 

% classsta is the observations vector O for the first trial of the word 
‘Statistics’. = 
classsta=vec2ind(ys); 
ys2=sim(net,csta2’); 
classsta2=vec2ind(ys2) ; 
ys3=sim(net,csta3’); 
classsta3=vec2ind(ys3); 
ys4=sim(net,csta4’); 
classstatest=vec2ind(ys4); 
ym=sim(net,cmicrol’); 
classmicro=vec2ind (ym) ; 
ym2=sim(net,cmicro2’) ; 
classmicro2=vec2ind(ym2) ; 
ym3=sim(net,cmicro3’); 
classmicro3=vec2ind(ym3) ; 
ym4=sim(net,cmicro4’) ; 
classmicrotest=vec2ind(ym4) ; 


yi=sim(net,cinstril’) ; 
classinstr=vec2ind (yi); 
vVi2—-Sim(net,Ginstr2“); 
classinstr2=vec2ind(yi2) ; 
Vius=sim (netyainstr3 4); 
classinstr3=vec2ind(yi3); 
yi4=sim(net,cinstrtest’) ; 
classinstrtest=vec2ind(yi4) ; 
vVoO-Sim(net ,cprotl’)-: 
classprof=vec2ind (yp); 
Woz=sim(net,cprof2 ’); 
classprof2=vec2ind(yp2) ; 
Vpo=sim(net,cprof3’ ); 
classprof3=vec2ind(yp3) ; 
yp4=sim(net,cproftest’); 
classproftest=vec2ind(yp4) ; 


% we divide by the constant multiplied before 
cstal=cstal/convert;csta2=csta2/convert;csta3=csta3/convert;csta4=csta4 
/convert;w=w/convert; 

enero !—-CMmicroly Convert ;, CmicroZ2—cm1icro2/ convert: cmicro3=cmicro3/convert 
Pemicro4=cmicro4/convert; 
Ginstril=cinstri/convert;cinstr2=cinstr2/convert;cinstr3=cinstr3/convert 
;Cinstrtest=cinstrtest/convert; 

SpLofl=cprotl/ convert; cprof2=cprof2/convert;cprof3=cprof3/convert;cprof 
test=cproftest/convert; 


% Test of vector quantization 
figure (1) 


uO prOt( 4, 2,0) plow: Gyecten (lee). *2 1-8 wiclasssta(1),:)),title([’ve 
ctorl,Class=’ num2str(classsta(1)) ’ M=’ num2str(M)])) 

Plubplou(4, 2,2) ;pon(l c,estal(2,-)>°*°,1:8;wiclasssta(2),:)),title([’ve 
ctor2,Class=’ num2str(classsta(2))]) 

BUDpHOt(4,2,35) plot (1: Grestalts,:),°*',1:8,w(classsta(3),:)),title([‘ve 
ctor3,Class=’' num2str(classsta(3))]) 

SUDDPIOE (2 72) plore Gl-8,cstal(4,:),°*’,1:8,w(classsta(4),:)),titie(|[‘ve 

) 


Gtor4,Class=’ numZstx(classsta(4))]) 
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Ssubploe(4, 2,5). plo (36, estate: 
ctor5,Class=’ num2str(classsta(5) 
SUbpIOE (4,2,6) ;plot (1:8, cstar(e,- 
ctor6,Class=’ num2str(classsta(6) 
Sup lot.(4:.2,.7 plot (ie 6vestal(). - 
ctor7,Class=’ numZstr(classsta(7) 
SHipp lot (4, 25-3) plot (ies nesta i( 3.,.: 
ctor8,Class=’ num2str(classsta(8) 
figure (2) 

SHipp lot (42,1) prow wioremi cro! (1,5) 7 
[‘’vector1,Class=’ num2str(classmicro(1)) ' 


eee ee eee eee ee Se 


— ~~ — 


J 


Sulbplot(4, 2,2); piouene ce cmicrol (2,7) ) —* 4: 


[‘vector2,Class=’ num2str(classmicro(2))]) 


subplot (4,2, 3) premier. So emicrol(3,:),°*. 2: 


[~vector3, Class-— suum str(classmicro(3)):))) 
Subp lot(4 2 4s promis -cmicrol (4,>),.° ==; 
[‘vector4,Class=’ num2str(classmicro (4) ) ] 
SUbDpLOE(4, 225) )5 lot (5, cmrerol (5,2 )9) 77. 
[‘vector5,Class=’ num2str(classmicro(5)) ] 
Silo pa oie 4eeztro.) Oo NOt Wi G, CMuerOl(6,.-) 40 * i 
[‘’vector6,Class=’ num2str(classmicro(6))]) 


Suplomucwe ys DlLOt(1:8,cmicrol(] «=: )) "2.1; 


[{vyeetoms, class=" “mumZstr(classmicro(7)})) ])) 


Snippets, 6) -ploc (+8 ,ecmicnoums,:),°** 7 1: 


[’vectors,Class=’ num2str(classmicro(8))]) 
sda(l:ia.)—wliclassmicro(i:T) ,<:); 
$qa(1:T,:)=w(classsta(1:T),:); 


xe oe 
Ao ron 
ee eo 


er oe 


w(classsta(Sje)) titled ve 


w(classsta(6),:)),title([’ve 
w(classsta(7),:)),title([’ve 


w(classsta(8),:)),title([’ve 


28,wi(classmicro(!),:)), rele 
M=’ num2str(M)]) 

G67w(classmicro(2) 7) )jerelet 
8,w(classmicro(3),:)) ,titlel 
:8,w(classmicro(4),:)),title( 
+o 7wielasSsmicro(S). 2), clemet 
16,wiclassmicro(6)7%)), ue re. 
8,w(classmicro(7),:)),title( 
8,w(classmicro(8),:)),title( 


bee en 
Gistances (1)=) Gist (cstalin, :) wiclasssta(n).) 
ayStancemi(nr= aise (cmicrol(n-:)) wiclassmicroin),:) )- 
end 
figure (3) 


stem(distances),title(’Euclidean Distance original/quantized vectors’) 


figure (4) 


stem(distancem) ,title(’Euclidean Distance original/quantized vectors’ ) 
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APPENDIX A2. FEATURES EXTRACTION/VECTOR QUANTIZATION USING 
-K-MEANS 


% Filename: trainingl.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Features Analysis (LPC & Energy) of 4 words (Statistics, 
MEerOosOne;s INSLCrUcEOr, 


% and Professor), 3 trials each + 1 test word, vector 
quantization using k-means 

% for HMM classification use. 

clear 

M4 — oo; 

aad t=2 Ob 


rand(’seed'’,1); 


load stat.mat;load micro.mat;load instrprof.mat; 

load words.mat; 

% signals interpolation 

%$ T is the number of segments 

% sstal,ssta2 ... coefficients of each trial 
[cstal,T]=coeff(stal) ;[csta2,T]=coeff(sta2);[csta3,T]=coeff(sta3);[csta 
4,T]=coeff(sta4); 

fentero ir) —Ccoekrt (microl) > [cmicro2, T]=coefrt (mienoZ); 
fcmuecros, 1 | —coechimimicro3); [cmicrod, Ff) —coefrt (micro4) ; 
fetmcesl, T)i=coeffi (instr!) [e1nstr2, Til=cockt (inserz): 
feamstrs, 1 | =coefi(anstrs) ; [Cinstrtest; 1 =ccetiimstrtest): 
feprorl  T|=Coectf(pEoEt): (cprotZaiijqecoctt(protZ):- 
[cprof3,T]=coeff(prof3); [cproftest,T] =coeff (proftest) ; 


%$ vector quantization 

pl=(cstal;csta2;csta3; 

SMlGrOl Gmicroz»CmiCros ;CinStrl- ecinstrZz;Cinstr si 
% CodeBook creation: 

[CODE Tabet, dist] = svqiplyeM, nett); 

% Vector quantization 

[w classsta, dist] = vq(cstal, CODE, n_it); 


iw2;-ClasssStaZz, Gise| = vo(cstaz; CODE es neiL); 

ifs Classstas. dist) = va(csta3, CODE, meic) > 

[wt, classstatest, dist] = vq(csta4, CODE, n_it); 
iw classmicro, dist] = valemicrol, CODE, M21); 
iveclassmicroz mat stilm=ssel(emicro2, CODE, nit); 

fer classmicroes, Gist) = Velemicros,, CODE, nit) ; 

iy -claSsmicrotest, Gist) ="volemicro4, CODE, n_it) ; 
fer eclassinstr,dust| = vel(emsrerl CODE, n_it) ; 
Iweeclassinstr2, aise) = voleinsterZ2, CODE, nm -it): 

lite Glassinstr3, G@ise| ="Vvalcinstr3, CODE, n_it)-; 

in; classinstrtest, Gist] = valcinstrtest, CODE, n_it); 


iw -classproft;, dist! = yvoleprotl,. CODE, n_it) : 
iw Classprorz atstl = VWaleprorc,, CODE, n_it); 
[weclassprors Gist le—) Vor cemors ) CODE, N_it) > 
[w classproftest, dist] = vaq(cproftest, CODE, n_it); 
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% Test of vector quantization 
Elgure |) 

Ssubolot (4, 2,1)splor( ls 7eeraimelia.) 
‘ num2str(classsta(1)) ‘° M=’ numZstr(M) )) 


Semloeloc(4,2,2);pletl(l- 8, escakt(2,:),°* , i: 
‘ num2str(classsta(2))]) 
Sumo lot (412.15) ;OUOEl ts wecrea (3, :).° rE 
‘ num2str(classsta(3))]) 
SUbplot(4, 2,4) plot eeestal(4,:) 7 1 
’ num2str(classsta(4))]) 
SNpplot (472, >) -plorwes cstal(57:), "724 
‘ num2str(classsta(5))]) 
Swboploe (4, 2,06 > lotr eo -cstal(6..)7°* ar 
‘ num2str(classsta(6))])) 
SHbDDLOE (47277) 2 plore: 6 -Cstal (7 a2) 7p? * 7) 
‘ num2str(classsta(7))])) 
sulee Pot (4723.6) pier! :8,cstal(8,:),°*",1 


‘’ num2str(classsta(8))]}) 

figure (2) 

Subp LOE) Ol Ot (1) 7 omemme sour ( Ii) 
[*vectorl Class=" numZstr(classmicro(])) 


rer ines 


val 


S,wtii,title([’vectorl,Class= 


8,w(2)),title([’vector2,Class= 


>8,w(3)),title([’vector3,Class= 
>8,w(4)),title({’vector4,Class= 
-3,wl5))yesele (vectors, Clasce 
:8,w(6)),title([’vector6,Class= 
:8,w(7)),title([’vector7,Class= 


:8,w(8)), title([’vector8,Class= 


rO7wiclassmicro()) =) )frxmelet 


7M emum2serm (My) 


Subptor 4a 2,2) -olot (1:8, cmrerol (2,3) * 2 ee awielassmicro(2Z) 71) ete len 
[*vector2, Class=' numZstr(classmicro(Z)))) 
Shipper as | 7 DLOL IL: 8)cmicrol(3,:),’ *" i isnwielassmicro(3) 72), eileen 
[‘vector3,Class=’‘ num2str(classmicro(3))]) 
Suppletii2, 4) -ploc( ls ,emiecrol (4,2), 47 ode cawiclassmicro( 4), 2) ) ee len 
[‘vector4,Class=’ num2str(classmicro(4))]) 
Sop kot 4.2, 5);pLOee ¢-cmicrol(s7:), -%"> bee wielacsmicro(5),2))7eieles 
[*“veetor5,Class=" numZ2str(classmicro(5)9)) 
Subploe (4) 2,6) > ploemes cmierol (6,2). = ls 8 7wielassmicro (6), 2) reutle 
[‘’vector6,Class=’ num2stri(classmicro(6))]) 
Slbplot (4,2, 7); ple@rmim@es ,cmilcrom? >)" = 1.0 wlelassmicro(7),:))perele 
[vector, ,Class=" mnumZstr(classmicro(7)) )]) 
Sibploel4, 2,0) -ploc( lo, cmacroelia,:), *°; les,wirelassmicro(8),:)),titlet 


[’vector8,Class=’ num2str(classmicro(8) )] 
tqatiet,:)J=w(classmicro(i: 7), :); 
$ga(1:T,:)=sw(classsta(1:T),:); 


fOr l= le 
distances (n)= -aisticstalin,:-)3wlelassstavn ,:)); 
aGistancemin)= dist (cmicrol (mn. ewicrleassmicrvoin), °)): 
end 
figure(3) 


stem(distances),title(’Euclidean Distance original/quantized vectors’) 


figure (4) 


stem(distancem) ,title(’Euclidean Distance original/quantized vectors’) 
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APPENDIX A3. FEATURES EXTRACTION (LPC + ENERGY) FUNCTION 


a 


% Filename: coeff.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: fuction that computes the LPC coefficients and the energy of 
each segment 

%$ T: total number of segments 

% ec: Coefficients vector 


function [c,T]=coeff (word) ; 

% Signals interpolation 

$word=word-mean (word) ; 

t=length (word) ; 

% interpolation of word, so every word has 10000 data points 
iMib-intewo. (1 -awora, @- iy 10000; 10000). 

% Pre-emphasis filtering: 

mre=-fB1lter([1,—.98],1,int); 


start=1; 

step=1200; 

1=10000-step; 

% Overlaping of segments: 

overlap=round( (10/100) *step); 

%# of Observations T: 

T=l1/step; 

T=floor (T) ; 

fOr nals Tt 
c(n,1:8)=ar_covar(pre(start:start+step-l+overlap) ,7); 
start=start+step; 
x=xcorr(pre(start:start+step-l+overlap) ); 
lx=length (x) ; 
center=find(x==max(x)); 
ec(n,1)=x(center) ; 

end 

Ce. y= een) 7 max (CeCe...) )= 

Qc, |) =e se -5- 
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APPENDIX A4. FORWARD VARIABLE ALGORITHM IMPLEMENTATION 


% Filename: forward.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Implements the forward algorithm (Rabiner 1986 pg9) 


function afw=forward(a,b,pi,0O,N,T) 
O: the observation sequence 
T: total number of segments 
a: the state transition probabilities (N x N) 
b: observation probability distribution (N x M) 
Di I dtetabestacemarstr1oution (N x 1) 
N: total number of states 
M: total number of possible observation symbols (dimension of 
odebook) 
afw: the forward variable (N x T) 
bbw the backward variable (N x T) 


GP dc 1 oP cP A dP oP oP oP 


% initial values: 
afw=zeros(T,N); 
ow ile oie: Ni) a td oN Oe 
% recursion: 
% we are using eps (for zero results) to avoid divided by zero problems 
if norm(afw(1,:))<eps 
afw(1,:)=eps; 
end 


% FW algorithm: 
for t=1:T-1 
Eom a=) 


S=0; 

foLrei—i 

% summation S: 

S=Stafw(t,1)*a(li,j); 
end 
$logafw=1lo0g10(S)+1lo0g10(b(j,0(t+1))); 


afw(t+1,j3)=S*b(j,0(t+1)); 
if norm(afw(t+1,j))<eps 
afw(t+1,j3)=eps; 
end 


tafw(t+1,j3)=10*logafw; 


end 
end 
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APPENDIX A5. BACKWARD VARIABLE ALGORITHM IMPLEMENTATION 


% Filename: backward.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Implements backward algorithm (Rabiner 1986 pg9) 


Puncemonmobw—-backward(a, 0,01, O0,N,2P) 


the observation sequence 

total number of segments 

the state transition probabilities (N x N) 

observation probability distribution (N x M) 

i initial sstate distribution (N scm) 

total number of states 

>: total number of possible observation symbols (dimension of 
codebook) 

% afw: the forward variable (N x T) 

% bbw the backward variable (N x T) 


dP dP oP dP oP OP OP 
= 2, 0 Of) © 


% BW Algorithm: 
% Initial Values: 


bbw=zeros(T,N) ; 
Dow(T, 1 N)=1-: 
% recursion: 
% we are using eps (for zero results) to avoid divided by zero problems 
1f norm(bbw(T,:))<eps 
bbw(T,:)=eps; 
end 


for t=T-1:-1:1 
for i=1:N 


S=0; 

Om waa oN 
S=S+a(1,3)*b(j,0(t+1) ) *bbw(t+1,3); 

end 

bbw(t,1)=S; 

1f S<eps 
bbhw(t,1)=eps; 

end 


end 
end 
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APPENDIX A6. FW/BW VARIABLES SCALING 


% Filename: scale.m 

%$ Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Implements scaling of FW and BW variables, according to 
Rabiner(1989 pg 272) 

% We perform scaling to avoid exceeding the precision range of the 
computer 

O: the observation sequence 

T: total number of segments 

a: the state transition probabilities (N x N) 

b: observation probability distribution (N x M) 

pi: initial state distribution (N x 1) 

N: total number of states 

M: total number of possible observation symbols (dimension of 
codebook) 

% afw: the forward variable (N x T) 

% bbw the backward variable (N x T) 


de dP dP GP dP dP GP 


function [afwl,bbwl,c]=scale(a,b,pi,0O); 


N=length(a); 

S=size(O) - 

T=s(2); 
afw=forward(a,b,pi,0O,N,T); 
bbw=backward(a,b,pi,0O,N,T) ; 
afw2=zeros(T,N); 

atwZtie p=alwtl,:): 


G(l)=l7sumtatw(l,:)); 
abwl (l,i 2N)=e(1)*afiw(i, i:N); 


%$ summation S: 
S=S+tatwl(t=1,3)*a(j,1); 
end 
aiwe (ce) 1)=S sot, Ot eb); 


end 
su=sum(afw2(t,:)); 
if su== 
su=eps; 
end 
c(t) =1/su; 
aewil(e,.)=C(t) “atw2 (ta 
end 
bbwl=zeros(T,N) ; 
bewin i. ).=e(1); 
fore c= v1 — 1: | 
tome a=. lh) 
S=0); 
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Fone a= 
S-Ss1ay) Sg, © (tt) bowl (e+1, 4) ; 
end il 
Dewi t= Siac te) - 
end 
end 


Sl 


APPENDIX A7. MULTIPLE OBSERVATION SEQUENCE 


% Filename: newmodel.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Builds the O_multi matrix, which contains all trials 
obsrvations and 

% calls trainmulti.m function 

%$ T: total number of segments 

a: the state transition probabilities (N x N) 

b: observation probability distribution (N x M) 

pi: initial state distribution (N x 1) 

N: total number of states 

M: total number of possible observation symbols (dimension of 
codebook) 


dP dP dP d® o® 


% models re-estimation 

% initial conditions: 
O_multi=[{classmicro;classmicro2;classmicro3]; 
$O_multi=[classsta;classsta2;classsta3] ; 
$O_multi=[classinstr;classinstr2;classinstr3]; 


%O_multi=[classprof;classprof2;classprof3]j; 


(ayo pie tralnmules (Oomulti,N, 1,M); 


APPENDIX A8. BAUM-WELCH ALGORITHM (MODEL REESTIMATION) 


— 


% Filename: trainmodel.m 
% Written by: M. Zambartas 
%$ Date Last Modified 10 August 1999 
% Purpose: this function implements the Baum-Welch reetimation 
algorithm as discribed 
in Rabiner (1986) pg 11 
the observation sequence 
total number of segments 
the state transition probabilities (N x N) 
observation probability distribution (N x M) 
1: initial state distribution (N x 1) 
total number of states 
: total number of possible observation symbols (dimension of 
codebook) 
% afw: the forward variable (N x T) 
% bbw the backward variable (N x T) 
% afwl: Scaled forward variable (N x T) 
% bbwl: Scaled backward variable (N x T) 


dP dP dP dP dP dP dP oP 
= a 0 0 MH] 


function [a,b,p1,P] = trainmodel (0O,N,T,M) 

% if the model is ergotic keep a=rand(N). If left-right then 
a=triu(a), because it has 

% to be upper triangular, in order to converge to an upper triangular 
matrix 


% Initial estimation of model lamda=(a,b,pi): 

$a=rand(N) ; 

b=rand(N,M); 

a=triu(a); 

% we need the summation or all rows of a and b to be 1, so: 
tm. D/.0(Ssum(b.)) .°*ones(i,M))> 


em ae, © (Ssumia.”’)).' *ones(1,N))- 
ae tame ee ty 

Pam =27 (50-4); 
end 


% we need the summation of pi to be 1: 
p= pi, /sum( pi) - 
for times=1:20 


eiw=Lrorward(a,b,pi,0,N,T); 
bbw=backward(a,b,pi,0O,N,T); 

%$ Scaling: (Rabiner 89 pg272) 
% afw scaling: 

anoLd=—a ; 

Eeold=b: 

Piola =p x - 


afw2=zeros(T,N); 
atw2 (i se=atwill.:)-; 


G(i)= 7 stm iarwil,:)): 
arwilte Ni) =c(l) *atw(i, 1 3N); 
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% summation S: 

S=Stafwl(t-1,3)*a(j,i); 
end 
atw2(t, 1) =S bli7O(e)).. 


end 
su=sum(atw2(t, 2): 
% to avoid ‘divided by zero’ errors: 
if su== 
su=eps; 
end 
e(t}=l/su- 
atwl (tb, j-eumieotw.(t, >) 
end 
% bbw scaling: 
bbwl=zeros(T,N); 
bowl) =) ed 1). 
£Or t—T—-1-—1:1 
foie oy 
S=0; 
fee ec) 
S=Stavieg)* bli, Obs!) ) *bowl(t+1,3)? 
end 
PS wieer ao eee 
end 
end 


% a reestimation 
in@ ae) 
S2=1; 
fOr ee — 1. — 1 
Sl=(.); 
if Oise d— 1) 
Si=(SGie arwl(t, |) *ali, 4) “bowintt ie) lg, O(t+1)))- 
end 
S2=(S2° sum(S1)]; 
end 
dena (i) =sum(sum(S2)); 
end 


fe Oise eed 
for j=l1:N 
SO oe 
£Or t=12T-L 
Sl=(Sl; afwllts1) “as bevel “bia Ott. |) i 


end 
numa(i,j)=sum(S1); 
1£ dena(i)== 
dena(i)=eps; 
end 
a(i,j)=numa(i,j)/dena(i); 
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wardt 


E 


03 


ies | 


fis 


3 
ee 


: oF. 
pay afi 


end 
end 


if norm(a-a_old)<.01 _ 


a=a_olda-, 
end 
% b reestimation 
fom 3 —ih 
Por c=—1 3M 
S$1i=0;S2=0; 
FOrSc=1 +E 
Si=Si -atwaree, a) *“Dbwitt, a) ~ (Ove) —=k) += 
S2=S2+afwl (t,j3) *bbwl(t,j); 
es — (0) 
S2=eps; 
end 
Digek)=Sl/s2- 
end 
end 
end 


if norm(b-b_old)<.01 


b=bEo la: 
end 


% P_bw(O|lamda) probability 

logP=-sum(log10(c)); 

P=10*1logP; 

pat p= (alwl ee). Dow (1 .2)./C( 1) P- 
ZlOgp1—(loglOtarwl (ly 7) )+logl0 (bbw (1,72) )—logil0 te(1) ) ) -loge, 
$p1i(1)=10*logpi; 

Lee WoOmmGol—D1 Old) =) UL 


onl=jonk ler 
end 


$fprintf(’%g %g %g \n’,norm(a-a_old) ,norm(b-b_old) ,norm(pi-pi_old) ) 
if norm( (a-a_old) ==0) &(morm(b-b_old) ==0) & (norm(pi-pi_old) ==0) 

break 

end 

end 


hail. 
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APPENDIX A9. BAUM-WELCH ALGORITHM (MODEL REESTIMATION) 
WITH MULTIPLE OBSERVATIONS 


% Filename: trainmulti.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: this function (re)estimates the model (a,b,lamda) for 
multiple Observation 

% sequence according to Rabiner (89) pg 273 

% Each row of O_multi is each Observation O, where: 

%* O: the observation sequence 

T: total number of segments 

a: the state transition probabilities (N x N) 

b: observation probability distribution (N x M) 

Di> initialstate distribution (N x 1) 

N: total number of states 

M: total number of possible observation symbols (dimension of 
codebook) 

%$ afw: the forward variable (N x T) 

% bbw the backward variable (N x T) 


dP oP oP HP HP oP 


Puncmtonml ao, bt = trainmulti(O multi Nt 
% s=size(O_multi); 
Erlals—s (i) - 
aEwiltmulen=([ |; bowl multi=|]-; 
% Computing P(O,lamda) & scaled fw and bw probabilities for each Obs 
(trial) 
fOr k=l triare 
[a7b; pl 72) = "trarnmode! (O_multr(k7:) Nii); 
Pamultitky=e- 
fatwihsbowlesel—scale(ia,b,pi,O multi(ky.) )- 
arwlemilci=larwlo mult. atrwl |; 
bowl muiler=(bpbwl multx bowl]; 
end 
% Model’s re-estimation using an alternative formula: 
$P_multi(l)=1; 
$a_old=a; 
ee epe st SIL SI) 
for j=l 3N 
Sik=() Ss2k-(); 
£Or "k= see rals 
Sl=1e s2= (4). 
for t=)-7—1 
S1=[(S1 afwl_multi(t,1i+(k- 
) *N) *a_old(i,j)*b(j,0 multi (k, t+1)) *bowlL multi (t+1,3+(k—1) *N) 1; 
S2=(S2 atwlemubemie, 14 (k-1) +N) *bbowl multe (t,1e(k-1)*N) le 
end 
Sl=sum(Sl)/7 P_mu Ler 
S2=sum(SZ)/P mulemGe. 
S1k=[(S1k S1]}; 
S2Zk=(SZk S2ae 
end 
S1k=sum(S1k) ;S2k=sum(S2k) ; 
a(i,j)=S1k/S2k; 
end 


oP oP dP oP cP cP cP OP CP CP F CP CP CP oP CP OP 


86 


$end 


% a denominator evaluation: 
fot lle Nl 
Sk=[]; 
for k= bra ls 
S2-—()5 
foe tale 
S1=(]; 
POmeg— —>N 
Sl=(Sl> atwilomulea(t aS e=1)-N) “ata = bbwl_ multi (eriga ek 
ies) je Olomu Leek, t+1))); 
end 
S2-—(S2>sum(Si)]; 
end 
Sk=/SkeS2/P_multi(k)]> 


end 
dena (i)=sum(Sk); 
end 


% a final Re-estimation 
Ore = 1: Ni 
for j=1:N 
Sk=[]; 
Poiamic= ida d als 
Si=([)5 
fer t=i. T=) 
Si=[S1 afwl_multi(t,i+(k-1) *N) *a(i,j) *bbwl_multi (t+1,j5+(k- 
ie pty ,Oemulei(k tei) ) )-: 


end 
SkaiSkesum (Si) / Pom ta (ik) |; 


end 
a(i,j)=sum(Sk)/dena(1i); 
end 
end 
% b reestimation: 
ioe i= 1 > N 
Oi eee 
Sik=[];S2k=[]; 
for k=l; trials 
S1i=[();S2=[]; 
LOY st=1 20 
Si=[S1 afwl_multi(t,j+(k-1) *N) *bbwl_multi(t,j+(k- 
1) *N) * (O_multi(k,t)==1) ]; 
S2=[(S2 afwi_multi(t,j+(k-1) *N) *bbwl_multi(t,j+(k-1) *N) J; 


end 
Si =Sum( Si) Pome o2—sumi S2)/ Pomulti.(k) - 
Sik=([Sik Si; S2k—(S2ZK.S2 |; 
end 

Sik=sum( Silk) SZ k=sum¢ S22); 

big ,1)=S1k/S2k; 

end 

end 
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APPENDIX A10. VITERBI ALGORITHM 


% Filename: viterbi.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: function that implements the Viterbi algorithm as discribed 
an 

% Rabiner (1986) pg 11 

the observation sequence of testing word 

total number of segments 

the state tranSition probabilities (N x N) 

observation probability distribution (N x M) 

1: Mnitialeystare wai1stribution (N x 1) 

total number of states 

: total number of possible observation symbols (dimension of 
codebook) 

% afw: the forward variable (N x T) 

% bbw the backward variable (N x T) 


oP dP dP dP dP dP oP 
s20OM HO 


fincCrEMOMmLEi,mi, 6 seq])=Viterbi(a,5,p1,0,N, 2) 
% Viterbi algorithm: 
tara valiie £11 (1)=pi(i)*bi (OL) 
fi=Zeros (1..N) >; 
Panella oat Lon CLs N,OCL) ) 5 
% backtracking pointer mi: 
Mme=Zeros (in N) ; 
Porm 2. 
fOr.) —1 2 N 
M=[]; 
fora cy 
Maer — yd) ati, 3) | 
end 
ArgMax=find(M==max(M)); 
fimeeme —masc (NM). ~b(7,0(L) )3 
Mime y= aroMax (1) ; 
end 
end 
% Path (state sequence) backtracking: 
s_seq=zeros(1,T); 
% final state @ T time: 
U=£1nd (fa Cie —— ee ee =) > 
s_seq(T)=u(1); 
for t=T-1:-1:1 
s_seq(t)=mi(t+1,s_seq(t+l1)); 
end 
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APPENDIX A111. SCORING (CLASSIFICATION) 


= 


Filename: score.m 
Written by: M. Zambartas 

Date Last Modified 10 August 1999 

Purpose: Evaluates P_bw(O|lamda) and P_viterbi(0|lamda) 
O: the obsServation sequence of tested word 

T: total number of segments 

a: the state transition probabilities (N x N) 

b: observation probability distribution (N x M) 

Di anitial state distribution (N x 1) 

N: total number of states 

M: total number of possible observation symbols (dimension of 
codebook) 
% afw: the forward variable (N x T) 
% bbw the backward variable (N x T) 


dP dP dP oP OP OP dP dP oP dP OP 


function [P_bw, P_v]=score(a,b,pi,0O) 
% Scoring of an Observation 
N=length(a); 

s=size(b); 

o— length (0); 

% Viterbi probability: 

f2t7mi, Sseseqi=viterbi(a,b,p1,0,N,1); 
Pevy=masc( f(t. )) : 


% Baum-Welch probability: 
abw=LOrward(a,b,p1,0,N,2) ; 
bbw=backward(a,b,pi,0O,N,T); 
{fafwl,bbwl,c]=scale(a,b,pi,O); 
% Rabiner 89 pg 273 P(o,lamda)=1/prod(c) 
sc=-sum(logl0(c)); 
P_bw=1/prod(c) ; 
if P_bw==NaN 

P_bw=0; 
end 
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APPENDIX A12. CLASSIFICATION RESULTS PLOTS 


Filename: scoretest.m 


Writt 
Date 
Purpo 


_bw(0| 


en by: M. Zambartas 

Last Modified 10 August 1999 

se: Plots classification results for all words, evaluates 
lamda) and P_viterbi(0O|1lamda) 


for every model 


Ox 
oes 
a: 
b: 
ae: 
Ne 
a 
odeboo 
afw: 
bbw 


the observation sequence of tested word 
total number of segments 
the state transition probabilities (N x N) 
observation probability distribution (N x M) 
initial state distribution (N x 1) 
total number of states 
total number of possible observation symbols (dimension of 
k) 
the forward variable (N x T) 
the backward variable (N x T) 


O=classstatest; 


‘ 


test: 


Statistics’ 


[P_bwl, P_v1]=score(a,b,pi,0O) 
f (P_bwl1==0) 
P_bwl=eps; 


i 


end 
Prey i=—()) 


i 


Pay 


end 


=eps; 


O=classmicrotest; 
‘test: microsoft’ 
[P_bw2,P_v2]=score(a,b,pi,0) 
f (P_bw2==0) 

P_bw2=eps; 


a 


end 
EP y2—-—0) 


a 


P_v2 


end 


=eps; 


O=classinstrtest; 


c 


test: 


Instructor 


[P_bw3, P_v3]=score(a,b,pi,0) 
EP pw3——0) 
P_bw3=eps; 


al 


end 
f- (Pav ——O) 


i 


S 


P_v3 
nd 


=eps; 


O=classproftest; 


f 


test: 


professor’ 


[P_bw4,P_v4]=score(a,b,pi,0) 
fee 0) 
P_bw4=eps; 


ct 


end 
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it (ea ——0) 
P_v4=eps; 
end 


subplot (2,1,1),, barh(10*log10([P_bwl P_bw2 P_bw3 P_bw4])); 
title([’Baum-Welch Probability - 1:Statistics, 2:Microsoft, 
3:Instructor 4:Professor, M=’ numZ2str(M) ’, N=’ numZ2str(N)])), 
xlabel(’‘’P_B_W[dB]’), ylabel(’tested word’) 


Craig 

SUbD MeL Zon) aebari O~LoglO{(Payvir Pp v2 P v3 P_v4])); 

peetele viterbl Probability — 2@2Statistics, 2:Microsoft, 3: 1lnsecuctor 
eeraressor, —M=  mumZser(M) “7 N=snumzstr(N)))), 
Rlabel(’P_V_i_t_e r b.3/dB)]"), ylaber@ tested word’ ) 

grid 
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APPENDIX A13. VECTOR QUANTIZATION ALGORITHM 
Ponceron [CODE, Label, distlw= svq(x, “ev ner): 


$svq Vector quantization using successive binary splitting steps. 
% Use: [CODE,label,dist] = svq(X,lev,n_it). 

% The final codebook dimension lev should be a power of two. dist 
% returns the distorsion values at the end of intermediate step. 
%$ n_it is the number of iterations performed in each step. 


% Version 1.3 
% Olivier Cappi, 28/09/94 - 04/03/97 
% ENST Dpt. Signal / CNRS URA 820, Paris 


% Needed functions 

if (exist(’vq’) ~= 2) 
error(“Funeetom vq 1S Missing.” ); 

end; 

% Turn verbose mode off 

QUIET = 1; 

% Input agruments 

@rroru(Margenk(3, 3, nargin)); 

% Dimension of imput data 

(nep)  ="size (xX); 

% Number of spliting steps 

Moese=—eaound (log2 (Lev) ); 

lev = 2*nbs; 

% Fixed perturbation 

Derturou— 0.01; 


% Initialize first centroid with global mean 
CODE = zeros(lev, p); 

CODE_ = zeros(lev, p); 

CODE (3...) = mean (xX) ; 

label = ones(n,1); 


fOr 1—i.nbs 
% 1. Codebook splitting 
fore j—i2 (ad )) 


CODE (2a), -) = (l+perturb) * CODE (3,2); 
CODE (24-775) = (l-pertume)) =. CODE(a,.)); 
end; 
%$ 2. K-means optimization 
Meme 2257) label, vdist) = vealx7CObe. (1:2°1,. 5) ,nZ21e, QULET):; 
Gist: (1 )ib=)Vaist (n_ it); 


BOiGgtnik le “COCeDOOK.Size %d:\te.sfim 42 °-1,01st (1) ); 
end; 


APPENDIX Al4. VECTOR QUANTIZATION USING LGB ALGORITHM 


function [CODE_n, label, dist] = vq(X, CODE, n_it, QUIET); 


vq Vector quantization using the K-means (or LBG) algorithm. 
% Use: [CODE_n,label,dist] = vq(X,CODE,n_it) 

% Performs n_it iterations of the K-means algorithm on X, using 

% CODE as initial codebook. 


% Version 1.3 
% Olivier Cappt, 28/09/94 - 16/07/97 
% ENST Dpt. Signal / CNRS URA 820, Paris 


error(nargchk(3, 4, nargin)); 
Hemtiarginy <.o) 

QUIET = 0; 
end; 


% Dimensions of X 

fae pie= Size(x)= 

% Codebook size 

m = length(CODE(:,1)); 

% Initialialize label array 
label = zeros (lyn) ; 

% As well as distortion values 
east (= zeros (1, n it); 


% Main loop 

CODE_n = CODE; 

mor Vtehe—aie-n 1b 
% 1. Find nearest neighbor for the squared distortion 
DIST = zeros(m,n); 


Gan (ome 8) 
eG) 1g ee = eee 
Dicotwiven. = SUMUL x — Ones i(n, |) *CODE MiG.) 2): 
end; 
else 
% Beware of sum when p = 1 (!) 
DIST =| (ones (ml) *x* =— CODE n*ones(1,n))2 72; 
end; 


{vm, label] = min(DIST); 
% Mean distortion 
dist(iter) = mean(vm); 
$2. Update the codebook 
meow = iO 
Owe, 1 em 
tc 9 = (dom) es 
ind = ind((label == 1i)); 
te (lengthi tne) )—— 0) 
% Isolated centroid are not modified 
MoOuk = hMeour = i 
elseif (length(ind) == 1) 
% When there is only one nearest neighbor for a given codebook 


entry 
CODE =n i: = ax ine 2. 
else 


25 


CODE snifey)) —-smean (|X Gama). 
end; 
end; 
% Affichage 
1f (~QUIET) 
SIeLoOrintt (1, Lteraktion sd Vet 38 \n 7 iter are e ier) 
end; 
Po e(neout > 0) 
% fprintf(1,’' Warning : %.0f isolated centroids\n’,n_out); 
end; 
end; 
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APPENDIX B. SEISMO-ACOUSTIC SONAR PROJECT INFORMATION 


This appendix contains some basic information regarding the NPS seismo- 


acoustic sonar project. Further details may be found in [10, 11, 12]. 


2D 


APPENDIX B1. SEISMIC WAVE ACTUATOR 





Actuator with waterproof case and coupling device [12] 
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APPENDIX B2. BEACH TEST SITE 





so 


Re 2 gal 
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Beach test site with data collection equipment [12] 
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APPENDIX B3. MINE-LIKE OBJECTS 





Powder keg target with lid open to show access [12] 





Gas cylinder target [12] 


98 


APPENDIX B4. BURIED MINE-LIKE OBJECTS 





Buried powder keg target with top removed [12] 


CORE NUE GEG G8, ste ag 


the, > gasiememner Petree tS 
Awine Thy ax? om, 





ie 
as ae ooo ee 
se es ate 


~ an Fra fy, 
a % , 


GOs ok, gs OR RA. MD 
. ed hates ot 


Buried cylinder with end cap removed [12] 
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APPENDIX B5. 8-CHANNEL DATA PLOT [12] 





Channel 1 Channel 2 
@ a 
CO SS A lop el rere re ee ee ace 
= = nt ; 
= Bqel----- 1__ File > tgt24lbsum 
& = 10Fites | 
@ @ 
TS an 
= = 
= = 
Ee = 
i fl 
@ @ 
TS a 
= = 
= = 
= c 
< mi 
@ a) 
Sey a=; 
= = 
oS cat 
= S 
= mi 
Range[it] Range[ft] 


8-Channel data plot of the received signals. Channels 1 and 2 come from the 
accelerator meters. Channels 3, 4, 5 represent the x, y, z motion of the 1% geophone, 
and channels 6, 7, 8 the x, y, z motion of the 2"° geophone. 
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APPENDIX B6. TEST SETUP FOR HELIUM GAS TANK AND GUNPOWDER 
KEG TARGET TESTS (INCREASING MASS) [12] 


This appendix contains the experimental setup for target detection with increasing 


mass tests conducted on November 6, 1998 and November 10, 1998. 


Source #1 
Orientation: Vertical 
Voltage: 

Frequency: 80Hz 
Cycles: ] 
Geophone #1 

Filter: High Pass 40Hz 
Gain: 40db 


10V (20Vpk-pk) 


Source #2 
Orientation: Vertical 
Voltage: 1OV (20Vpk-pk) 
Frequency: 80Hz 
Cycles: ] 
Geophone #2 
Filter: High Pass 40Hz 
Gain: 40db 


Comments: Overcast, med-high tide, 4-5ft waves. 
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Helium Tank (Nov 6") 


Gunpowder Keg (Nov 10") ‘ 


APPENDIX B7. Rayleigh wave [11] 


S-Wave R-Wave 











Ww 


rE (+ down) 







Particle Motion 
u 


(c) 


Direction of Wave Propagation 
es 


Seismic wavetrain resulting from a single vertical impulse source [11] . 


APPENDIX B8. GEOPHONE [11] 





Triaxial geophone seismic sensor [11]. 
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APPENDIX B9. CROSS POWER [12] 


Geo #2 Real Power 


Geo #1 Real Power --> tgt4lbsum 


ee re ¢ @ @ @ 


1 
t 
t 
1 
rf 
’ 
’ 
’ 
’ 
7 
i 
’ 
1 





apnyyd uy 


Imaginary Power 


Imaginary Power 


1 
1 
1 
+ 
1 
1 
1 
1 


1 
\ 
! 
' 
+ 
! 
\ 
! 
! 
T 
! 
! 
t 
Lf) 
O 





apnyydury 





apny|duy 


Range [ft] 


Range [ft] 


Real and imaginary cross power components of the received signal for both 


geophones. Target located at 22 feet 
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APPENDIX B10. TARGETS STRENTH 


Target Strength Vs Target Mass 


Target Strength in dB 


Total Cylinder Mass in Ibs 





Target strength vs. target mass for cylinder target [12] 


Target Strength Vs Target Mass 


ce 
"~~ 
= 
= 
SD 
= 
@ 
Pe 
= 
A 
= 
@ 
GD 
he 
cs 
Se 


Total Powder Keg Mass in Ibs 





Target strength vs. target mass for powder keg target [12]. 
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APPENDIX B10. IMAGINARY CROSS POWER SIGNAL FOR THE CYLINDER 
TARGET 


Vector Polarization Filter (Geophone #1) 


150lbs I 
254lbs 


358lbs 





moplitude 
© 


A 
© 
Q 


150lbs 
254lbs 
358lbs 
462\bs 
566lbs 


Amplitude 


Range [ft} 


Imaginary power plot for cylinder, November 6‘" experiment, 5 weight types [12]. 
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APPENDIX B11. IMAGINARY CROSS POWER SIGNAL FOR THE POWDER 
KEG TARGET 


Vector Polarization Filter (Geophone #1) 


120Ibs 
224|lbs 
328lbs 
432|bs 


o36lbs 


Amplitude 





21 22 23 24 25 26 2d 28 2g 30 
pale ay 
Vector Polarization Filter (Geophone #2) 


120Ilbs 
0.02 224lbs 
a 328lbs 

432|bs fm 


© 


536lbs 


Amplitude 


Caen) 
Oo CO 
te NI 





Range [ft] 


Imaginary power plot for powder keg, November 10" experiment, six data plots. 
Note that the 2" geophone did not receive the target [12]. 
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APPENDIX C. MATLAB CODE;-HMM BASED CLASSIFIER FOR MINE 
RECOGNITION 


This Appendix contains the various MATLAB files used for recognizing the 


mine-like objects using HMMs. MATLAB files written in [11] used for preprocessing 


the data are also included for completeness. 
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APPENDIX Cl. COMPUTATION OF THE CROSS POWER SIGNAL FROM 
RAW SIGNAL; 8 CHANNEL; SIGNAL. CODE WRITTEN BY M. 


FITZPATRICK [12] 
% Name: Present3.m 
% Author: LT Mike Fitzpatrick 
% Updated: 9/31/98 


% Description: This program conducts a Hilbert analysis.... 


$$é%$éInput Parameters%%% 


Range=1; $Turns-on plotting with range axis 
wavespd=295; *wavespeed [ft/s] 

t. stare=0 050; %Set start time [s] 

to stop=0-. la0F %Set stop time [s] 

Eo stare io- $Set start range [ft] 

r_stop=28; Set stop range [ft] 

scale=0.0230; Set axis scaling (Set to 0 turns-off scaling) 
georange=10; Enter range to geophone [ft] 


%%%Date, Directory, & File%t%% 
date=’Nov_6’; 
directory=’Target3’; 

cd (date),cd (directory) 


COUNT=0; $Set start count 
while (1) 
Eley eicp( ~““Hilbert s2nalysis. SUbrOoULIne =. Jair *.mat-; 
COUNT=COUNT+1 
1f COUNT==1,load tgt0Olbsum, end 
tf (COUNL=—-2, Loac, taoutibsum, end 
ie COUNL—=—=3)loqdGg tots osum-. end 
if COUNT==4,load tgtlZ2lbsum, end 
1£ COUNT==5, load tgtl6ibsum, end 
if COUNT==6, break,end 
[M,N])=size(channel) ; 
transform=hilbert (channel) ; 
Pwr (2 2 COUNTL) Cony (transtorm( =, 37555 erales -orm 2,5): 
Pwr (22, COUNT) —con) (tCranstorm( +67) -transtrorm(:, $8); 
end 
end 


%%%Set Axest%% 

if Range==1 
[maxl, index1]=max(abs (channel (:,1))) 
[max2,index2]=max(abs(channel(:,2))); 
index=round ( (indexl+index2) /2); 
start=round( ((indexl+index2) /2)+(r_start/ (dt*wavespd) )); 
stop=round (((indexl+index2) /2)+(r_stop/ (dt*wavespd) )); 
range=wavespd* (t(start:stop)-t (index) ); 
$start=round(delay+(r_start/ (dt*wavespd) )); 
$stop=round (delay+ (r_stop/ (dt *wavespd) )); 
$range=wavespd* (t(start:stop) -t(delay) ); 

else 
start=round( (t_start/dt)+1); 
stop=round((t_stop/dt) +1); 


e 


/ 
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end 

1£ scale== =i 
Realmax=1.1*max(max(abs (real (Pwr(start:stop,:))))); 
Imagmax=1.1*max(max (abs (imag (Pwr(start:stop,:))))) 

else 
Realmax=scale; Imagmax=scale; 

end 


bon u—1 :COUNT—1 
Pwrl(:,n)=imag(Pwr(start:stop,1,n)); 
Pwr2(:,n)=imag(Pwr(start:stop,2,n)); 
end 


SEBPLlOCCings%% 
figure,orient portrait 
if Range== 
Subp Loe (2711 71) 
pillomlrange, Pwrl(:,.), eG , ‘LbineWidth’,1),hold@ 
plot (range, Pwrl(:,2),‘m’, ‘Linewidth’ ,2 
ELOEt(range, pwr :,3) ,°q°, “Linewidth 3 
plot (range, ewril(:,4),’b’, ’LineWidth 74 
ploE( range, Pwrl(:7.5)> Yr’, ’Linewidth’ 7s) 
axis([{min(range) max(range) -Imagmax Imagmax]) ,hold 
title(’Vector Polarization Filter (Geophone #1) ’) 
xlabel(’Range [ft]’),ylabel(’Amplitude’),grid 
legend( t56lbs’, 160lbs’ , 364lbs’ , 4268ilbs‘*, 57Z21bs’) 
EFLEEEEEEEEBESES 
subplot(Zyl, 2) 
PlLOt (range ,Pwr2(:,1), °c’ >’ binewidth’ -1) 
plot (range, Pwr2(:,2),’m‘’, ‘Linewidth’ , 2) 
3 31), 
4) 


‘ 


~ 


) 
) 
) 
) 


f 


;nold 


PlOet(randerrwre(:,3)5 9°, Lanewidth’ 7 
plot(range, PwrZ(= 4), bb’, LineWreen’ , 
plot (range, PwrZ(:,5), 4r©°, Linewidth’ 75) 
axis([min(range) max(range) -Imagmax Imagmax] ),hold 
title(’Vector Polarization Filter (Geophone #2) ’) 
xlabel(’Range [ft]’),ylabel(’Amplitude’),grid 
legend( 156lbs’, ~160lbs’ ,’3641bs’ , “4684bs’ , S72bs’ } 
end 


jel. 
oc. . 


Lil 


APPENDIX C2. DATA SELECTION USED FOR THE MULTIPLE WEIGHT 
SIGNAL EXPERIMENT SET-UP. DATA USED FOR HMM TRAINING 


% Filename: kegcyl.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

%$ Purpose: Generates and plots the signals for the multiple weight HMMs 


training 

% kegmat.mat: imaginary cross power of the keg signal 

% cylmat.mat: imaginary cross power of the cylinder signal 
% back.mat: imaginary cross power of anon target signal 


load kegmat 

load cylmat 

load back; 

range=linspace(1,60,3000) ; 

% interpolation to 3000 points of all vectors, in order all signals to 
have the same 

% length 


Keg4 32 Mowery) =) | > kegqasolb (2277...) = Iepeyl> seb (2277, )= 1h 
sk=size (keg2241b) ;sc=size(cy13641b) ;sb=size(back) ; 

kegint=zeros (3000, 6) ;cyllint=zeros (3000,6) ;backint=zeros (3000, 6) ; 
% 6 trials per signal: 

£or n=1:6 


kegzzatibine(:0)=(interpl(1:sk(1), keg2241b(- om) ie sk(1)/300l=ski 2); 
keg43 7 oine(s,n)=(interpl (1:sk(1)), keg43215(:,n),1-sk(1) /3001:sk(l)));: 


kKegsS Olbimt (29n)=(Gnterp! (lesk( lj kegs olb(: mn) -sk(1)/3001l-sk()) a; 


cylse4tbine(-,) = (intere |) (lwcek e673 64 lbw) sc(l)/300l: seth); 
CyVl46S 1 Dime (2 ie (Imcerpl( 1] cei) cyl etoomet. mm) 1.sc(1)/3001:scay yy 
eviS7 2 bint (Gem (antcerpL( Pseq vey los7Zio-7 nn). 1-sc(l)/s001l:sc{ lj) 


backint (+0) —(imterpllh-coo Gl) back ae sb(l)/3001-sb(l) je] 
end 


% plot of all signals 
SUDDEGE (S274) 
plot (range, keg22 4 bint ec | GInewidth’,1),hold 
plot{(range, keg2241bint G2) om . Ginewidth’,1) 
plotl(range,Keqg224Ubine Geo 6 oo banewiath” , 1) 
plot (range, keq2 24 "bine 4 >.) eeameWwr ath, 71) 
plot (range ,keg2241bink (se kh Fe eanenidten’, 1) 
plot (range, keg2241 bint. oe. =. eimewidth”’, 1) 
hold 
title(’Powder Keg 224l1b ‘'’) 


axis lOv6o0. -0.037 0.05) \2oqmid- 
SUbpmor (S62, 3) 


Blot range, Kegeoel bint (.) l) .7e">) Linewidth’,1)),hold 
plot (range, keq4321bant(:, 2) 7am |)  LineWidth’ , 1) 

plot (range, keg4321lbint(:,3),'@, ’LineWidth’,1) 

plot (range, keg4321lbint(:,4),’b’, ‘LineWidth’ ,1) 

Dlot (range, keqosc bint (:,5)5°k 3 Linewidth’ 71) 

plot (range, keq432lbint(:,6),’r’, “LineWidth’ ,1) 
axis({0 60 -0.03 0.03])),grid; 

title(’Powder Keg 4321lb ’) 


ie ke! 

SubpLottse. 5) 
plot (range, keg536lbint(:,1),’c’,’LineWidth’,1),hold 
Plot( range, kegssolbint(:,2);,'m ~ binewideth’ ,.1) 
plot (range, keq536lbint(:,3),’9’, Linewidth’ , 1) 
plot (range, kKeGoJtolbint(:,4),"b | Linewrden 71) 
plot (range, keg536lbimem@,5),’°k’, LinewWidth’ ,1) 
DIlOtirange Kegosololnt.:., 6), LL --binewlatn ; 1) 
esos | 060 --0-05 0-03] grid; 
title(’Powder Keg 536lb ’) 


hold 
SubplOor(s.2, 2) 
DIGE (wange, Gyso2lbint (+, 1),'c’s Linewidtnee. 1 )ehold 
DlOottrance, ©vyrso4ibint(:, 2), ’' mm’, Linewidtm.1) 
pict (range yey L3641bint(:,3),°9q",’ LineWidth 71) 
plot (rangeyey] 364 lbint(:,4),’b’, ‘LineWidth’ ,1) 
DIOE(Lrande7cy 3564 bint(., 5), kK’) Linewidth 2) 
Bloptrange,cy.so4l1bint(:,6), ©, Linewidrn 71) 
eecus tO 'o0 —-0. 02 -0.02));qr1d) 
Piltlet .cviander 3641b *) 
hold 
SubpLetwee a4.) 
plot (sange,cy M4Geciternt (+71), ‘ec, LineWidth 71),nhold 
Dlot(mange, cylaoslbinet.,2)), mm’) Lbinewidth:; 1) 
DLOM Grange, cyl oclbintl:) 3). 9 > Linewidcth |.) 
Dice (range, Gy leoclbint(:,4)) bb, Linewidth 77) 
DIOL (range, cyl4oetbint(2,5)7°k’, “Lbinewidth’ 1) 
Dl@E(Langer cy 466 lbint(:76)9 m7 Linewidth-, 1) 
exis ( (0 60 -020270.02))), grid: 
title(’Cylinder 468lb ’) 
hold 
Subplor(4s, 275) 
plot (range, cy57 Zaieant(:,1), cc’, “Linemadth’,1),hold 
plot (range, c~yloS/Zibint(:,2),'m’, ‘LineWidth’ ,1) 
DloOti(range, GYVls72tbine(-73),°C', Linewiaen”,1) 
plot (range,cy 5 /Zibint (2,4), b' ." Linewrden’ , 1) 
DlLOe (Tange, cy 572) bink( 2, 5). K's binewidtn’ , 1) 
DlGti range, cvls?2hoime (7,6). °F’, “manewidth’ ,1) 
reas ct © 160-0 .02 0.02)), grid; 
Peele cylinder S721b-* ) 
hold 
%$ target location: 
for m—l:6 
1k224 (n) =find 
(keg2241bint (1000:1800,n) ==max (keg2241bint (1000:1800,n))); 
1k432 (n) =find 
(keg432 Wpame( H000-1200>n ) ==max (keg432 li bint (1000-1800,n))): 


NS 


e536 (i) =find 

(kego36lbint (1000:1800,m)—=maxtkeq536lbint(0G0- le0Gr im). 
1k224 (n) =1k224 (n)+1000-100; 
1k432 (n) =1k432 (n)+1000—106- 
1k536 (n) =1k536(n)+1000-100; 

end 

mee 2Z4=71k224( 1) 54 k432=uka ee) i ko 36=i1koSo) 

Bor n=1:6 
ies 6 (hh). & ere) 

(cvl364 ibint ( LOOG ser n)==max (cyl3641bine VO00 sO na 
1¢364 (n) =1¢364 G1) +1000-100- 
1c468 (n) =find 

(cyl468 1 bint (1O00F Fe00) n) =—max(ev1 468 loink (100 EeOtr rns 
i1c468 (n) =ic468 (n) +1000-100; 

1eo0 2 ae 

(evi57Ztbpint (000. Ve007n)=—max (evil 5721 bint (VO00 eGo ms ae. 
1ea7 Zin) =165724n) 1000-100; 


end 
1C304=16364 (1) -1c468=1C468 (1) -16572=16572 (6) 


% so now we know where the target is located, we can build the signal 

file 

% signal (data,mine, trial) 

% mine(1-3) is for the keg (6,8,10 ft)’ s target, mine(4) is for the 

cylinder’ s target 

% and mine(5-10) background and other non target data 

Signalmulti=zeros(400,10,6); 

for trials=1:6 
signalmulti(:,1,trials) =keg224lbint (1k224 :i1k224+399, trials): 
signalmulti(:,2,trials) =keg432] bint (1k432 :1k432+399,trials); 
Signalmulti(:,3,trials) =keg536lbint (1k536:1k536+399,trials) ; 
signalmulti(:,4,trials) =cyl1364lbint (1c364:1¢c364+399,trials); 

Signalmudte1( 3, 5,Erials) =ev 14681 bank ve468— tet6e1+399 Exials)- 

Signalmulti(:,6,trials)=ecyloe7 ZI bint ites) 22165724599; frials): 


Signalmulti(:,7,trials) =backint (16364 2163645599; trials): 
srqnalmilti( -7e,trlalsS) =pbackint(tk536- ess o75 Io erials) - 
sigmanmulti(:,9 trials) =backinte (2001-24007 tridails); 
signalmulti(:,10,trials) =backint (1500:1899,trials) ; 
Signalmulti(:,11,trials) =backint (800:800+399,trials); 
signalmultit(:, 127trials}) =baekine (2500-2699, erials)- 


end 

$signalmulei (: 37) =Signalmuliim= 5,2) 72; 

for m=1:4 

for n=1:6 
Signalmulti(:,m,n)=stenmalmiltit: m,n) mean (Ssignalmilts (mon) 

end 

end 

ea 

save signalmulti signalmulti 

cd data 

figure (2) 

hange=1:400; 

$range=linspace(21,29,400) ; 
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SuUoerouto, 2, .) 
plot (range, stonalimulta (: laa ec 7 Linewidtn’ 71), hold 
plot (range, signalmulti(:,1,2);/“m’, ‘LineWidth’,1) 
DLOu (range, > onalmuleit: les) qo tanewidth’, 1) 
DlLoOL (range, > nemalmuLed.(:,1,4) 5) bo Linewiden’ , 1} 
DlOt | rance crmamaimuier (+,1,5),°k*) “Linewidth’ , 1) 
pDloul(range-stane multi: ,1),6)oe 2, Linewidth’ ,1} 
title(’Powder Keg 224lb '’) 
hold 

SupDpLOb (5,270) 
Bum range, stonalmidti(:, 2,1) eee Linewidth’ ,1),hold 
PLoeywrange, signalmulti(:,2,2), Mm) Linewidth’ , 1) 
ploetrange, signalmulta(:,2,3), Ga". bainewidtn’ ,1) 
pmOe (range s1ggamultci(:,2,42),'> —} binewidtn’ , 1) 
ploe(range, signalmulti(s,2,5),° kk", LineWidtn’ ,1) 
plot (range, signalmult ime, 2,6), ’°x', Linewidth’,1) 
title(’Powder Keg 4321lb ’) 
Hoke 

Som lot (3; 2;.5) 
Bloc (ranoge,>+reoMalmuLtit:, 5,1); °c", Banewidth’ 71),hold 
DlOt (range, stddaimuilta (3,2), ’ m ;, Lainewidth’ ; 1) 
DlOtlrance, o1ddalmulei(:73,3), 6 , Linewiden 1) 
DlOttrance, slonalmiler(: 75,4). bb’) binewidth: . 1) 
DloOttrancge,sionalmulei(.,3,5)5. k’, Linewisdthn- >) 
DloOtteancge, sltonalmulti(: 3,6), ©, baneWwidtm- 4) 
title(’Powder Keg 5361b ’) 

axis({0 400 -0.02 0.02}) 
xlabel (’Data Points’) 

no ld 

SubplOt (3,2, 2) 
pigttrance,stonalmulti(:,4,1)4 "e", “Linewidth’, 1) , hold 
pugtiuange,si1gonalmulti(;:,4,2),’'m’, 'Linewidtny’ , 1} 
mrot (range,sSigqgnalmulti(:,4,3),%9q’ , ’Linewidth’ ,1) 
plot (range, signalmultei<:,4,4),°b’ ,’ Linewidth’, 1) 
PrOE(vange,signalmulti(: 4,5), Kk’, ‘Linewidth’ ; 1) 
DLOt (tange, See atm 6) Go einen arnaa 
title( “Cylinder 3641b *) 
nella 
SuUpplOwls, 2,4) 
DlOt (mange -etonalmultit 5, 1),°e", “Lanewiden’, 1), hold 
pDloul(range, stgnalmulti(: 5,2), mm’ 5 * binewidteh’ , 1) 
plot (range, signalmulti(:,5,3),’g’, ‘Linewidth’,1) 
DlOE (mange, sitegmea limiter: 5,4) .°D" ,“lamewidth” ;1) 
plot (range, stonalmulti(:,5,5),°k’,’ Lbinewidth’ ,1) 
MUO (range, stonalmulel(:,5,6), 6, LineWidth’ ,1) 


Elile( Cylinder 4681p =) 


hold 
Seoplot (3,2,6) 

plot (range, signalmulti(: os hee fe lbainewrdth’1),hold 
plot (range, signalmulti(: Zw o, > Linewiathn.,1) 

plot (range, signalmulti(: on So. Lanewilath.— 1) 
DlOuVnanicge, sagnetmaltit =) 6,4), “b’, ’“LineWidth’,1) 

DUO (sancgerstonmalmudel(:760,5),'k’, “LinewWidth’ ,1) 
DlOe(ranGgessitonaimilti(:,6,6),’r’,’Linewidth’ ,1) 


cicle(’Cyiinderes7Z2lb ) 


les) 


xlabel (’Data Points’) 
hold 
ed. 2 


save signal signalmulti 
ayel elehe 
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APPENDIX C3. FEATURE ANALYSIS FOR MINE LIKE OBJECT SIGNALS; 
MULTIPLE WEIGHTS EXPERIMENT SET-UP 


% Filename: trainingkegcyl.m 
% Written by: M. Zambartas 
% Date Last Modified 10 August 1999 


% Purpose: Features Analysis for mine-like object signals with 
multiple weights 


% (3 sets of weights for the keg and the cylinder - 6 trials 
each, and 7 background signals), for 
% HMM recognition use. We will create two classes, one for 


every mine-like object. 


% Vector quantization using k-means 
% M: # of symbols 

% N: # of states 

%$ mines: # total # of mines 

SC: # coefficients matrix for all mines 
% w: # VQ coeff matrix 


% vector quantization with multiple lb (keg,cylinder, background) 
clear 

load signal 

signal=signalmulti; 

clear signalmulti; 

M=8;N=4; 

rand(’seed’,1); 

ei — 2.0% 

$segs=4; 

mines=12; 


s=size(signal);sl=s(1); 


%seg=fix(sl/segs); 
Coefficients evaluation 
$c=zeros(T,4,8,segs,6); 
for mine=1:mines 

Bor trial=1:6 


[e(:,2?,Mine, trial), ©) =coefeimterpisignal(:,mine,trial)); 


end 
end 


% vector quantization 
pl=[]; 
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for mine=1:mines 
foOreerial=)] 26 


Dil=[pl -cl s/s mine, cera 
end 
end 
[CODE, label, dist] = svea(ol, M, n_it); 


class=zeros(T,mines,5); 
w=zeros(T,8,mines,5) ; 
for mine=1:mines 
EOre thle h: 6 
{(w Class (7mane, trial), disse) =vaq(cel.,: mine, erial) , Coben iti 
end 
end 


% Test of vector quantization 


$figure (1) 
subplot 427) olet(1:8,c(1,:,1,1),** 2 8 wlelass(1 17 eee 
eC[”*vectorl, Class=" numZstriclass(171,1)) ~ M=" Vnum2str(M))) 


Ssoubplot(4 252));ol00(1:8,6cliwew 2). et erwielascs ly 12) 2 2 eaten 
e([’vector2,Class=’ num2str(class(1,1,1)) °’ M="* num2str(™) ]) 
Saqa(i: 7 ,.:)=wtelassmicro(1]:T),:); 
$qa(1:T,:)=w(classsta(1:T),:); 
distances=[]; 
$for n=1:segs 
Set oOummer ali; 5 
% distances= [distances 
Grsttcw-ese nn, criall wl(Cclass(: nn, Erial),:7n, trial) ) |; 
$end 


$end 
%figure(3) 
$stem(distances),title(’Euclidean Distance original/quantized vectors’) 
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APPENDIX C4. SIGNAL DECIMATION 


% Filename: coeffinterp.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: fuction that performs a decimation of segment and normalizes 
ae, 

%* T: total number of segments 

%* oc: Coefficients vector 


PinecELoMmic, {i =COCLTEIMeerp (word) - 
worad=word-mean (word) ; 


$int=interpl(1:1,word,1:1/10000:10000) ; 


tpre=filter({1,-.98],1,word) ; 
pre=word; 

% hanning filter: 
han=hanning (3); 

$pre=conv (han,pre) ; 
l=length(word) ; 


$start=1; 

step=1/2; 

Sl =Z2U000- 

% 0% overlap 
overlap=round( (0/100) *step) ; 

%# of Observations T: 
T=1/step+(l*overlap/step%’2) ; 
G=tloor (if): 

Lor n=1:T 
c(n,1:10)=interpl(1:step, pre(1l+step* (n-1) -overlap* (n-1) :step*n- 
overlap* (n-1)),1:step/10:step) ; 

end 

Bets 1) =e eeman(e¢:.1).) 5 
ees; 1) =cl., L}- 5- 

% normalization: 

for n=lL:T 

e(n) 2) -cims. )-y max tabs (e(n.:))) | 

end 
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APPENDIX C5. HMM TRAINING AND SCORING FOR MULTIPLE 
WEIGHTS MINE LIKE SIGNAL DATA 


% Filename: sequence.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Performs all the sequence of HMM training and scoring the 
signals from multiple lbs mines 

% rotating every time the tested signalkegmat.mat: imaginary 
cross power of the keg signal 


% cylmat.mat: imaginary cross power of the cylinder signal 
% back.mat: imaginary cross power of anon target signal 


% testing sequence 
seq=[2 3 yo 
Hl 
al 
al 
a 
aie, 
for ctest=1';6; 
ts=seq(test,:); 
% newmodelkegcyl trains the HMM using the seq sequence of Signals 
newmodelkegcyl 
% scretestkegcyl scores the HMM for every testing signals 
scoretestkegcyl 
P_bkeg224 (test) =P_bwk224; 
P_vkeg224 (test) =P_vk224; 
P_bkeg432 (test) =P_bwk432; 
P_vkeg432 (test) =P_vk432; 
P_bkeg536 (test) =P_bwk536; 
P_vkeg536 (test) =P_vk536; 
P_bcy1364 (test) =P_bwc364; 
P_vcy1364 (test) =P_vc364; 
P_bcy1468 (test) =P_bwc468; 
P_vcyl468 (test) =P_vc468; 
P_bcy1572 (test) =P_bwc572; 
P_vcy1572 (test) =P_vc572; 
P_bback (test) =P_bwb; 
P_vback (test) =P_vb; 


NON NY W 
WWW BP A 
PPO 
U1 OV OV OV OF 


lee 


end 


APPENDIX C6. GENERATION OF HMM OBSERVATION SEQUENCE FOR 
MULTIPLE WEIGHTS MINE LIKE SIGNAL DATA 


% Filename: newmodelkegcyl.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Creates the multiple observation matrix O_multi for multiple 
observations 

% HMM training of multi lbs testing 

% ts: traininf sequence, assigned at sequence.m file 

% training sequence: 


%® keg (224,432,5361b) : 


ZOniser—(elassi 1, 6S (1))*'>class(:) 1) tsi2))*-class(:,1,ts(3))"-elass(: 
mises (4) Class Gam, tse >Class(22,ts(1))*;-;class(:,2,ts(2))" -class( 
RePes( sje Classi 2 ts (4)) selass(2)2,ts(5)) seclass(:,3,ts(2))-class 
eo es cee eGlass(: 5 ,tSi(3))" -elass (@:75,-ts(4)) *-elass(:,3,ts8(5)) we 
pee vlander (364,466, 5721b) 
Semilet—(classt:,4,ts(1))’-class(:)4,ts(2))”’;class(:,4,ts(3))~-class(:, 
ees ine Classe £5 (5) ) °sClass(:;, S7es(l)) -class(r,5tsi(2))*:class( 
Ao ES ts) Rete. Sees) Velass( 2, Ses (5) Je class -y76, ts(1))  -elass( 
MO, ES (2) ie elass(+26,ts(3))°"sCclass(- 46, cs (4) -class (976, tsi(5)) 1. 


feeb plea be ratnmule (OumuLti, Nera M):: 
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APPENDIX C7. SCORING FOR THE MULTIPLE WEIGHT MINE LIKE 
SIGNAL SET-UP 


%* Filename:scorekegcyl.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Scoring of all testing mine-like signals-multiple lbs case 


‘test: keg2241lb’ 


O=class(:,1,test); 
[P_bwk224, P_vk224]=score(a,b,pi,0’) 
tf (P_bwk224==0) 
P_bwk224=eps; 
end 
Li FOP2VK224=—0) 
P_vk224=eps; 
end 
‘test: keg4321b’ 
O=class(:,2,test); 
[P_bwk432, P_vk432]=score(a,b,p1,0’) 
Peer mows 2——0} 
P_bwk432=eps; 
end 
if (P_vk432==0) 
P_vk432=eps; 
end 
‘test: keg5361b’ 
O=class(:,3,test); 
[(P_bwk536, P_vk536]=score(a,b,pi1,0’) 
LE CPs pwikossb6——0) 
P_bwk536=eps; 
end 
PE | (Pevico36——0) 
P_vk536=eps; 
end 


‘cest: cylinder3641lb’ 
O=class(:,4,test); 
[P_bwc364, P_vc364])=score(a,b,pi,0O’) 
if) (P -pbwe364=—0) 
P_bwc364=eps; 
end 
if (P_vc364==0) 
P_vc364=eps; 
end 


‘test: cylinder4681b’ 
O=class(:,5,test); 
[P_bwc468, P_vc468] =score(a,b,pi,0’) 
f£ (P_bwc468==0) 

P_bwc468=eps; 
end 
PEO VCooo==0) 


P_vc468=eps; 
end 


“test -scylinders 7215 
O=class(:,6,test); 
[P_bwc572, P_vc572]=score(a,b,pi,0’) 
pie  bweoe2=—0) 
P_bwc572=eps; 
end 
Pi eave os 7 2—-— 0 } 
P_vc572=eps; 
end 


‘test: background’ 
O=class(:,11,test); 
[P_bwb, P_vb] =score(a,b,pi,0O’) 
Pee owo=-—0) 
P_bwb=eps; 
end 
if (P_vb==0) 
P_vb=eps; 
end 


figure 

SUbDIOL (6,2, test *2-1), barh(1l0*1l6ogl0((P_bwk224 Plbwk432 P obwk536 
P_bwc364 P_bwc468 P_bwc572 P_bwb])); 

Pe cest—— 

title([’ Pr[O|lamda(cylinder)] - 1:keg_2_2_4 1_b, 


num2str(M) ’, N=’ num2str(N) ]) 
end 
axis([-50 01 7]); 
ie Lest=—6 
end 
1f test==6 

xlabel(’P_B_ W[dB] ’) 
end 
if test== 
ylabel(’‘tested segment’ ) 
end 


if test== 
ylabel(’tested segment’) 
end 


isi. 

hold 

subplot (6,2,test*2), barh(10*log10([P_vk224 P_vk432 P_vk536 P_vc364 
P_vc468 P_vc572 P_wb])); 

$title([‘Viterbi Probability - 1:keg_224lb, 2:keg_4321b, 3:keg_536lb, 
feacylingersoO41 bo Gy linder 46615, 6:cylinder_572lb, 7:background, M=’ 
muMmests (MM) ) N=" nmumZstn(N) |), xlabel(’P_V_it_er_b_i[dB]’), 
ylabel(’tested segment’) 
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axis({=50 0 1 7)); 
Grad 
Lf test== 
end 
Lieetest== 
ylabel (‘tested segment’ ) 
end 
hold 
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APPENDIX C8. DATA SELECTION FOR THE MULTIPLE DISTANCE SIGNAL 
EXPERIMENT SETUP. DATA USED FOR HMM TRAINING 


% Filename: multift.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Generates and plots the signals for the multiple FT HMMs 
training 

% kegmat.mat: imaginary cross power of the keg signal 

% cylmat.mat: imaginary cross power of the cylinder signal 

% back.mat: imaginary cross power of a non target signal 


load multiftkeg 

load cyl6noe 

load back; 

range=linspace(1,60,3000); 

ReqoGgeizz/O. =) weolULei 22/76, -)=([1e 

sk=size(keg6ft) ;sc=size(cyl6noe) ;sb=size(back) ; 

kegint=zeros (3000, 6) ;cyllint=zeros (3000, 6) ;backint=zeros (3000,6); 


% interpolation to 3000 points 

fon i= leo 
Kegqé6itinti:5n)=(interpl (l:sk(1) 7; KeGgOLt( tym), 1 ssi) /S00Lssk(1)))’; 
KegG Seimei im) = (interpol lock) keqgttt: ume esky) 7s 001: sk(l)))”; 
Kequgrtineesn —(anterpl(lsek(l) 7 keg lL0rhG. mie ch) 7 S001 skit) ))’; 


CV UGneel ie) —(interpl(l-sell)7cylonoe(=,n) lece (1) /300l2se(1))e; 
backwinet:. 9) —(interp! (1 =sb(1) 7 back(: jmigee: sbi) /3001:sbil)))’; 
end 


ebploe (4, 172) 
plot (range keg6 teint (+, 4), CC’ , bineWidth 71) 4 held 
pilot (range. <eqaecint(:.2) .-m , LaineWwidth’,1) 
plot (eamcge, keqotctint(:,3),’9°,°’Linewadth’ oy 
DilOtl range rcedo Betmintc. 4). 6 , Linewider , 1) 
( tomo). 4, mineWwidtn _ 1) 
( (ero )e  “o Linewrarn 1) 


plot(range,kegoéftint 
plot(range,keg6ftint 
hold 

title(’Powder Keg 6ft (total:24ft) ’) 


mes (10 60 =0.-02 10.02 por: 

BmMibplot (4, 1)2) 
plot (rangerkegag teint (2). C¢ , Linewidtib; 1),hold 
DLOE (range, Kegoteint (32). mM’, binewidth’7 1) 
plot | range, <equmtant(+,5),. 0 — eanmeWwidth ,1) 
DlOe( range, Keqamemnlt (4) b , MAnew lath gly) 
DlOt (range, Kegs Peint(:,5),’k’ », Linewidth 72) 
Dlot(range, Kegs fein (-,6) 07" binenwidte’, 1) 
axis( (0.60 —-0 10260202) ,Ggrre- 
title(’Powder Keg 8ft (total:22ft)’) 


nold 
SUDO LOule nl. 35) 


mlot (range, keqlOEtInEl< 71) Cc’, * Linewiden 7) hold 
plot (range,kegl0ftint(:,2),’m’, ‘LineWidth’,1) 

plot (range, keqlOftine(: ,3)3o ,° LInewiden ..) 

plot (range, kegl0ftint (294), °’b’,’LineWidth’, 1) 

plot (range, kegqlOfLtinti >). kk’ ,’ Linewiden- ») 

plot (range, keqlOftame(- 76), rr’, Linewiden. 
axis(i0 60 -0.0290,02]). grid; 

title(’Powder Keg 10ft (total:20ft)’) 


hold 

range=1:3000; 

SubplOt (4 
plot (ranges ey lGnoeermte(:,1),’c’ >’ LineWwidth 72) hold 
plot (range, cy lGnoeint(:,2),’m’,*Linewidth 71) 
prot (wange,e7laneeint(:,3), 9°," Linewidtn’, |) 
plot (range, eylenoeint(:,4),’b’,’LinewWwidth’ 7) 
plot (range ve, lonoeint(:,5),°K’ ,*Linewidth’ 2) 
DlOt | wamge, ey lonOcint(.,6),'°r’ 7 Lbanewidth’ ~) 

axis(1053006--0-02 0.02)]) ,orid: 
Citle<teyiinder 1O0fE (total-=20ft)~) 

Hold 

for n=1:6 
ik6(m)—-ftind (keq6f£tint(1000:1800,n) —=max(keq6titint (1000-18007 n))); 
ico) -Himae | kKeqottinte (O00: 13800 7m) ——-max keqettint (1000-18007 n i.) 
1k10(n)=find 

(ceglGbeine (2000-1800, n)==max(keql0Etint (1000; 1800 .n))) ; 
TK6m) =1k6(n) +2000-100: 
1k8 (n) =1k8 (n)+1000-100; 
1k10 (n) =1k1i0 (n) +1000-100; 

end 

16 seo ies aes CL) tel OSLO): 

for n=1:6 
he Loin) =i ne 

(cy bénoeint (1000-21800,n)==max(cyl6énoeint (1000:1800,n) )) ; 
re LOtm) —1elO(m)+lo00-100- 

end 

ic10=1c10(1)-100; 

% so now we know where the target is located, we can build the signal 

file 

% signal (data,mine,trial) 

% mine(1-3) is for the keg (6,8,10 ft)’ s target, mine(4) is for the 

cylinder’ s target 

% and mine(5-10) background and other non target data 

Signalmulti=zeros (400,10, 6); 

for trials=1:6 
Signalmulti(:,1,trials) =keg6ftint (1k6:1k6+399,trials); 
signalmulti(:,2,trials) =keg8f£tint (1k8:1k8+399,trials); 
Signalmulti(:,3,trials) =kegl0ftint (1k10:1k10+399,trials) ; 
signalmulti(:,4,trials) =eylonecine (1610516104399 7 erials)-: 
Signalmulti(:,5,trials) =backint (i1c10:1¢c10+399,trials) ; 
Signalmulti(:,6,trials) =backint (1k10:1k10+399,trials) ; 
Signalmulti(:,7,trials) =backime (2001: 24007 crials)- 
signalmulti(:,8,trials) —backine (1500-4299, Crurals): 
Ssignalmultit:, 9, trials) =pbackime (600-200 -59S triacs). 
Signualmulei(:,10,trials) =backint (2500-2899 erials)- 


‘ 


end 


Sigmalmu kei see asa gnalmulti( : ,3¢ge 3: 

for m=1:4 

for yn=1-6 _— 
signalmulti(:,m,n)=signalmulti(:,m,n)-mean(signalmulti(:,m,n)); 

end 

end 

edi. 

save signalmulti signalmulti 

cd data 

figure (2) 

Beam e= aaa) ; 

Supeiet (4, 15/1) 


Plu@t (range, signa lmumbtate et 1) coc’ linewidth”, 1),hold 
uot (range, signa lmMiuLei( 2 eee me. lnmewLoth’ , 1) 
DlLet (range, signature 1.5 )> Oo ,’ Linewidth’ ,1) 
Plot (range, sronalmulti(:, 1,4), b’, laneWidth’ ,1) 
PUetlrangie, sional mMibei(o1) Siac) Binewidth’, 1) 
Plot (aange, si gnalmulti(;, 1 {6jer , Bamewidteh’ ,1) 
hold 

Subplot (471) 2) 
ProE(range, Signmalmmlti(:;2;i))"e¢ >" binewidth’ ,1), hole 
Dloel(sange,Ssiqnalmuilti(:,2,2)- mm, UIinewicun. )) 
pilec(range,stona lnulem(s, 2,3)" 6 ) banewi den’ ,1) 
Pe@eet(rance, stgnalmultiz(:,2,4), b’, lamewidth’ , 1) 
Plkeeirange, s1gnalmulen(2,2,5), Kk . bamewidtn’ ,1) 
plot (range, signalmulti(:,2,6),’r’,’LineWidth’,1) 
hold 

subplot (4,1,3) 
DlOEtrange, Ss 1 Gnalmultats, 5,1) Cc. Lemewiaden >.1) held 
DOE (range, Slonea musta > 73,2 )e. mM’) bemewadth’ ; 1) 
DLOElronce, saGmhalmuilud( 25,5). 6°, lWameWwidth’ , 1) 
DLOE Geange, slonalimublbit2,73,4)0 >” Linewidem 1) 
pIlOE(range, signalmubei(:,3,5)).°% . banmewaodth’ , 1) 
PlOet-ange, siOnalinutG1 (2,356). 2. banewaden 5 1) 
Ine 

Suoplot (4, 1,4) 
plot (range, signalmulti(:,4, Pace Lanewidth’, Pineda 


= 


) 

ean. Lamewideth’ , 
),’g’, ‘LineWidth’, 
) 
) 
) 


AL 

plot (vance stena lmumes(:; 4, i 
Al 

7 oe GAnewicatn 7i 
alk 

1 


4 
4 
plot (range,signalmulti(:,4, 
4 
4 
4 


~ 


plot (range, signalmulti(:, 
plot (range,signalmulti(:, 
plot (range, signalmulti(:, 
lnieylites| 


ieee TinewLlatn | 
 Laaewactia: 


~~ 


a 


f 


a) 


APPENDIX C9. FEATURE ANALYSIS FOR MINE LIKE OBJECT SIGNALS; 
MULTIPLE DISTANCES EXPERIMENT SET-UP 


% Filename: trainingkmmulti.m 
%$ Written by: M. Zambartas 
% Date Last Modified 10 August 1999 


% Purpose: Features Analysis for mine-like object signals with 
multiple ft 


% (3 sets of ft for the keg and one for the cylinder - 6 
trials each, and 6 trials of background signals), for 

% HMM recognition use. We will create two classes, one for 
every mine-like object. 

% Vector quantization using k-means 


$M: # of symbols 
% N: # of states 


% mines: # total # of mines 

ee # coefficients matrix for all mines 
% Ww: # VO coeff matrix 

clear 


load signalmulti 
Signal=signalmulti; 
clear signalmulti; 
M=c  N=a. 
rand(’seed’,1); 
ae — 20), 

$segs=4; 

mines=10; 


s=size(signal);sl=s(1); 


$seg=f1ix(sl/segs) ; 


$Coefficients evaluation 
$c=zeros(T,4,8,segs, 6); 
for mine=l1:mines 

for trial=1:6 


(ol hereto. Ti=coctLinterp (signal(: -mine,trial))-> 


end 
end 


% vector quantization 
pitt, 
for mine=1:mines 
for trial=1:6 
Dla lenet =): -mine, brral) is 
end 
end 


(GODEnelabeiedi st). = svai(pl,. M, samt) - 


class=zeros(T,mines,5); = 
w=zeros(T,8,mines,5) ; 
for mine=l1:mines 
hor ceria l= tc 
ibveclass( “mime, trial) sdlstl=vaql(c(:, *,mine, trial), CODE, n.1t): 
end 
end 


% Test of vector quantization 


$figure(1) 

Tetomuet a cee) pmom so, C(l, +) aie (ot eo wiclass(1,1,1);°,27 2s titl 
e({’vector1,Class=’ Be ici eee ’ M=’ numZ@str(M M)}) 
wsrooloc(47272) ploe(l-S cil, = 72) | o,wiclass(1,1,2),: 2 ete aae 
e([{’vector2,Class=’ num2str(class(1,1,1)) ’ M=’ num2str(M)]) 

Sda(l:T, -)=wtlelassmicro(1-T),:); 

$gqa(1:T,:)=w(classsta(1:T),:); 


distances=[]; 
$for n=l:segs 
Soeeror trial=i> 
% distances= [distances 
Gustci ny trial) ~wteclass(=,n, trial)72,n, trial) ”)]-> 
$end 


end 
%figure (3) 
stem(distances) ,title(’Euclidean Distance original/quantized vectors’ ) 
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APPENDIX C10. HMM TRAINING AND SCORING FOR MULTIPLE 
DISTANCE SIGNAL MINE LIKE SIGNAL DATA 


% Filename: sequencemulti.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Performs all the sequence of HMM training and scoring the 
Signals from multiple ft mines 


% rotating every time the tested signalkegmat.mat: imaginary 
cross power of the keg signal 

% cylmat.mat: imaginary cross power of the cylinder signal 

% back.mat: imaginary cross power of a non target signal 


% testing sequence for multi target distance 
% rotation of the testing Signal: 
seq=[2 3 5 
se 
iL 
i 
1 
eZ 
bOmecest—!-6; 
ts=seq(test,:); 


bo dN ND W 
WWW PP HP 
Ci Ov Oy OF OV ON 


SHUI 
et 


newmodelmulti 
scoretestmulti 

P_bkeg6ft (test) =P_bwk6; 
P_vkeg6ft (test) =P_vk6; 
P_bkeg8ft (test) =P_bwk8; 
P_vkeg8ft (test) =P_vk8; 
P_bkegl1l0ft (test) =P_bwk10; 
P_vkegl0ft (test) =P_vk10; 
P_beyl1l0ft (test) =P_bwc; 

P veyll0ft (test) =P_vc; 


P_bback (test) =P_bwb; 
P_vback(test) =P_vb; 


end 
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APPENDIX C11. GENERATION OF HMM OBSERVATION SEQUENCE FOR 
MULTIPLE DISTANCES MINE LIKE SIGNAL DATA 


% Filename: newmodelmulti.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Creates the multiple observation matrix O_multi for multiple 
observation 

% HMM training of multi ft testing 

% ts: traininf sequence, assigned at sequence.m file 

% training sequence: 


% initial conditions: 

% class (coefficient(1-T) ,segment (1-4) ,mine,mass) 

% keg (6,8,10£t): 

Oenthei—liclass(- 1 ts il) ) -class(?) lo tst2))  -class(:,1,ts(3)) sclasstie 
Ienecan) (Classi ts(5))) -eclass(*¢2,tsiL. -class(:,2,ts(2))"sclass(- 
mres(3 i -class(=,2,ts(4)) “:;class(:,2,ts(5)) -elass(-73,ts(1))’ selass( 
Moree (2) class ( 30, tsi(3))’ -class (2,350s(4)0 > class( 37 ts(5)) 41> 

% cylinder (single 10 ft) 

FOeMUle r= |(Class{(:,4,ts(1))*’;>class(:747Cs(2))-elass(:,4,tsi3)>" -class(: 
Pores (a0) -elass(:44,ts(5))*); 


fa, oom — ce carmmu lta (Osmulta Nr MM); 


ley 


APPENDIX C12. TESTING DATA; SCORING FOR THE MULTIPLE 
DISTANCES MINE LIKE SIGNAL SET-UP 


% Filename: scoretestmulti.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Scoring of all testing mine-like signals-multiple ft case 


‘test: keg6éft’ 
O=class(;, 1 ,tese) 
[P_bwk6, P_vk6]=score(a,b,pi,0’) 
if e  pwko—-0) 

P_bwk6=eps; 
end 
if (PP vk6=s0) 

P_vk6=eps; 
end 
‘test: keg8ft’ 
O=class(:,2,test); 
[P_bwk8, P_vk8]=score(a,b,pi,0O’) 
i heGeeowieo—-—0) 

P_bwk8=eps; 
end 
Lin Prvks==0) 

P_vk8=eps; 
end 
‘test: keglOft’ 
O=class(:,3,test); 
[PmowitOwe=viclOl=score(a,b,p1,07) 
if SCP Sb )0==0) 

P_bwk10=eps; 


end 

Lf (Pevkl0—=0) 
P_vk1l0=eps; 

end 


‘test: cylinder’ 
O=class (274, ese) 
(Pi bweeP ve |=score(a.o, pi 0") 
fre E bwe==U) 
P_bwc=eps; 
end 
ineeeve——O) 
P_vc=eps; 
end 


‘test: background’ 
O=class(:,8,test); 
[P_bwb; P_vb] =score(a,b,p1,0") 
i£ (P_bwb==0) 
P_bwb=eps; 
end 
Peete vo=—O) 
P_vb=eps; 
end 


Stlbplor(G, 2, eestes 2-4), barh(l0*logle([P-bwks P_bwkS Pubwkl0 P_bwe 


P_bwb])); 
he test—— 
fie Le (ie Pr(O|lamda_k_e_g) - 1:keg_6ft, 


2 KeGgeore 5 KeGeiw Ort, 4:cylinder,~cest : background, M=’ numZstr(M) °, 
N=snumastem(N), ). T=2 *)) 


end 
if test== 

xlabel(’P_B W[dB]’) 
end 
if test== 
ylabel(’tested signal’) 
end 


axis([-50 01 5)); 

spaniel 

subplot(6,2,test*2), barh(10*log1l0([P_vk6 P_vk8 P_vk10 P_vc P_wb])); 
$title([’Viterbi Probability - 1:keg_6ft, 2:keg_8ft, 3:keg_10ft, 


a cylLingden, cest: background, M=" numZ2str(M) ">; N= numZzstr(N) ]) 
if test== 

xlabel(’P_V_i_t_e_r_b i[dB]} ’) 
end 


fre test==3 
ylabel(’tested signal’) 
end 


exis({-50 0 1 5)]); 


grid 
hedd 
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APPENDIX D. MATLAB CODE; NEURAL NETWORK BASED CLASSIFIER 
FOR MINE RECOGNITION 


This Appendix contains MATLAB files used to recognize mine like objects using 


a Supervised backpropagation feedforward neural network. 
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APPENDIX D1. NN TRAINING AND TESTING FOR MULTIPLE WEIGHTS 
SIGNAL MINE LIKE SIGNAL DATA 


% Filename: nnlbs.m 

% Written by: M. Zambartas 

% Date Last Modified 10 August 1999 

% Purpose: Performs mine like object signal classification - multiple 
weights set-up 


% using a supervised backpropagation feedforward neural 
network 
S es; testing sequence used for the testing/training rotation 


% O_multi_keg: multiple observation of keg signals 
% O_multi_cyl: multiple observation of cylinder signals 


% De input vectors 
S30: target vector; 1 for keg, 2 for cylinder 
5 Sco output matrix 


SCG=|(Z203545o Gls 4.5 6-1°2 4 5 Gel 2. 35 Oe 4 6 28 aoe 
sc=[]J; 


EOretLest—l36 
ts=seq(test,:); 


Ounulimmereg—-i|class(:,]1,ts(1) -elass(: 1, tsit2Z))) 7 cllass(:,1,ts(3)0" -elas 
See tee ce Lass (- 5, ts (5)) ;class(:- 2,05 (135) -class(:,2,ts(2) cla 
SS (wees s ye selass (> 72,tsS(4))*selass (:72 > ts(5)) ~class(:,3;ts(l)} scr 
ASSO ote)) -class(:73,0S(3))” -class(-,3,ts (4) 17 :class(:,3,ts(5)) 18 


Ommpletecy|—l(Glass(:-~47ts (lj) classi! -4 ts(2))] class(:,4,ts(3)) setae 
Sl roeate) iy elass (F465 (5) class( soe es iyo ellass(:75,ts (2)07 7eta 
SSi-wiomestisy) class (+7 5-65 (4). class (-y or rous mw oeiclass(:;6, 6S (1) } een 
asst 6S (2)) -elass (:76,65(3))  -eclasst= oees (4) 2; class(:;6,ts(5) 12 
% generation of input vectors (matrix) p: 

p=(O multi. keg Olmulti_cyl]: 


Ca) be eel eee a eee ee 2 2 2 ce 22 2 el, 
% # of hidden layers:60. # of output layers:1 Functions used: logsing, 
purelin for hidden 

% and output layer relatively 

Met=newlt (tl 8: Lesi io0vll i wlegsrce = curelin }) treinim ): 
figure (1) 

net.performFcn=’sse’; 

net.trainParam.min_grad=le-20 

net.trainParam.goal=0.01; 

net.trainParam. show=10; 

net.trainParam.epochs=300; 

InerE-er)|]=train(net,p, £)- 

Ol=class(:,1,test) ; 

keg224=sim(net,0O1) 

O2=class(:,2,test) ; 

keg432=sim(net,0O2) 
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©3=-class(:,3Ftest) ; 
keg536=sim(net,0O3) 
O4=class(:,4,test); a 
cy1364=sim(net, 04) 
O5=class(:,5,test); 
cyl1468=sim(net,05) 
O6=class(:,6,test); 
ey 15 /2=Ssim(net, O65) 
O7=class(:,7,test) ; 
back=sim(net,0O7) 
se=(se; 
Keg ver keoqds7, kegs36, cvl1364,>cvi4698, cyl5/72, back]; 
figure (2) 
BoUbDlOr(?,., L6St);>, Stem. 7, sc), axistiiiees) 2.5)]) hold 
end 
EeNN Cutput plot: 
SuloplOe (3535, lL) pp LOe tl 6-=se (7.1) ), axis (fio 04 )) erele 


({‘’keg2241lbs’]),grid 
Stipmlomr(s, os. 2)ye DIOL 6) Se(+,2)), axis (i670 4) citle 
({’keg4321bs’]),grid 
SUS lOct Ss, oS) iploc (il: o., set: 3) )>vaxis{ [P6700 4), title 
eiekeas36lbs”’ |) ,qrid 


SUoop lott 3,3 74) ~solot (1:6, scl: ,4)), -axis( [1 soeOeo he title 
Mi evi2odlbs* |] jearid 

Bul lotts,3,5) @aplot (1s6, se(:,5)), axis({(([1 6. 0-4)]), title 
Cl cyl46slbds’ ] j¥erid 

Stop lor 3.3 Oe plot (deo, Sels6))> axis{([1l 6 O84 pe-catle 
ae GVio7 215s >|) orid 

SulpLobis, 5; /), pict (126, se(.,/))-. axis( ib 6 0 4) )yeacle 
ol’ background’ ]).,grid 
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APPENDIX D2. NN TRAINING AND TESTING FOR MULTIPLE DISTANCES 
SIGNAL MINE LIKE SIGNAL DATA 


% Filename: nnft.m 

% Written by: M. Zambartas 

%$ Date Last Modified 10 August 1999 

% Purpose: Performs mine like object signal classification - multiple 
distances set-up 


% using a supervised backpropagation feedforward neural 
network 
% ts: testing sequence used for the testing/training rotation 


% O_multi_keg: multiple observation of keg signals 


% O_multi_cyl: multiple observation of cylinder signals 
5 Dp input vectors 
Sout target vector; 1 for keg, 2 for cylinder 
Se Se. SUlEpUL Matrix 


Se0=(2,394.5.6-1583 4 5 6-1 2 4.5 6:1 2 Seem 2 3 4-6-1273 ae 
sc=[]; 

for test=1:6; 

ts=seq(test,:); 


Ormudervekeq=(class(;.1,tstl))”,class(:,1,cst2)) -class(:,17ts(3))7 elas 

Stee cs (4) seeselass(:.1,ts(5))>class(:,2,tsi(l))- ;class(:)2,€s(2)5 cla 

SS 2,6S(3)) sclass(:;2Z, eta popeGi sa) “eclass|:,/3,08 Gee oct 

ass(:,3, nee ciines( €S(3)0. class (.-37Es(4y —-class(<25, ts (Span 

0 Bnulbiimeey i =(class (:,4,ts(1))*selass(:74,ts (2) -eracss(:,4,t8(3) elas 
eG: 4 cst 4) Se a: CS Se ee 


% generation of input vectors (matrix) p: 
=(O%mltiekeg Oemulti._cy] |; 


| pigeemem meee 1 lt ee el IM 2 2 22 ele 
% # of hidden layers:60. # of output layers:1 Functions used: logsing, 
purelin for hidden 
%$ and output layer relatively 
Net=Newru lilies oleeo i POULIN logsig, « purelim jy “trainim’)); 
ragure (il) 
net.performFcn=’sse’; 
net.trainParam.min_grad=le-20 
net.trainParam.goal=0.001; 
net.trainParam.show=10; 
net.trainParam.epochs=300; 
fmet, er) =train(net, po, c)-¢ 
Ol=cllass(:,1,test)- 
keg6ft=sim(net,O1) 
O2Z-clascs(:,2, test); 
keg8f£t=sim(net,0O2) 
O3=class(:,3,test); 
kegl0ft=sim(net,0O3) 
O4=class(:,4,test); 
cyl110ft=sim(net,04) 
O5—ceClass(:,5, test) ; 
back=sim(net,0O5) 
sc=[sc; 

keg6ft, keg8ft, keglOft, cyl1l0ft, back]; 
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figure (2) 


end 
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